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ABSTRACT 

This paper introduces and describes the data cubes from GHIGLS, deep Green Bank Telescope (GBT) 
surveys of the 21-cm line emission of H I in 37 targeted fields at intermediate Galactic latitude. The 
GHIGLS fields together cover over 1000 deg^ at 9.'55 spatial resolution. The H I spectra have an 
effective velocity resolution about 1.0 kms“^ and cover at least —450 < ulsr < +250 kms“^, 
extending to ulsr < +450 kms“^ for most fields. As illustrated with various visualizations of the H I 
data cubes, GHIGLS highlights that even at intermediate Galactic latitude the interstellar medium 
is very complex. Spatial structure of the H I is quantified through power spectra of maps of the 
integrated line emission or column density, Nm. For our featured representative field, centered on 
the North Ecliptic Pole, the scaling exponents in power-law representations of the power spectra 
of A^hi maps for low, intermediate, and high velocity gas components (LVG, IVC, and HVC) are 
—2.86 ± 0.04, —2.69 ± 0.04, and —2.59 ± 0.07, respectively. After Gaussian decomposition of the line 
profiles, A^hi maps were also made corresponding to the narrow-line and broad-line components in the 
LVG range; for the narrow-line map the exponent is —1.9 ± 0.1, reflecting more small scale structure 
in the cold neutral medium (CNM). There is evidence that filamentary structure in the H I GNM 
is oriented parallel to the Galactic magnetic field. The power spectrum analysis also offers insight 
into the various contributions to uncertainty in the data, yielding values close to those obtained using 
diagnostics developed in our earlier independent analysis. The effect of 21-cm line opacity on the 
GHIGLS Vhi maps is estimated. Comparisons of the GBT data in a few of the GHIGLS fields with 
data from the EBHIS and GASS surveys explore potential issues in data reduction and calibration and 
reveal good agreement. The high quality of the GHIGLS data enables a variety of studies in directions 
of low Galactic column density, as already demonstrated by the Planck collaboration. Eully-reduced 
GHIGLS H I data cubes and other data products are available at www.cita.utoronto.ca/GHIGLS. 
Keywords: ISM: clouds - ISM: structure - radio lines: ISM 


1. INTRODUCTION 

The 21-cm emission line of H I is the most commonly 
used tracer of the three-dimensional structure of the dif¬ 
fuse interstellar medium (ISM). Original studies were fo¬ 
cused on structures and kinematics within the Galactic 
plane of the Milky Way (Burton 1976; Kulkarni et al. 
1982; Stil et al. 2006), and have been expanded to stud¬ 
ies of the vertical structure of the disk and, more gen¬ 
erally, gas at intermediate Galactic latitudes extending 
into the Galactic halo (Heiles 1976; Dickey & Lockman 
1990; Kalberla & Kerp 2009; Putman et al. 2012) 

The discoveries and insights gained in these studies 
have benefited tremendously from the all-sky LAB survey 
(Kalberla et al. 2005), which is beam sampled at 36' an¬ 
gular resolution. Stimulated by these results, the range 
of motivations for further H I surveys is quite sweeping. 
We have used the 100-m Robert G. Byrd Green Bank 
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Telescope (GBT, Prestage et al. (2009)) at the National 
Radio Astronomy Observatory (NRAO^) for complemen¬ 
tary studies at intermediate Galactic latitude, where of¬ 
ten the line of sight column density is low so that the 
emission is not strong. This paper is to introduce, de¬ 
scribe, and disseminate the data from our deep GBT 
surveys of the 21-cm line emission of H I in 37 fields 
at intermediate Galactic latitude, obtained mainly with 
the Auto-Gorrelation Spectrometer (ACS) over the pe¬ 
riod 2005 to 2010 (see Blagrave et al. 2010 for a prelim¬ 
inary report). We refer to the project by the acronym 
GHIGLS (GBT H I Intermediate Galactic Latitude Sur¬ 
vey). The total area mapped is over 1000 deg^ and al¬ 
though this comprises only about 2.5% of the sky, the 
judicious choice of environments sampled means that a 
broad range of scientific questions can be addressed, as 
outlined below. Compared to LAB, the GHIGLS data 
have a higher angular resolution of about 9' and are 
Nyquist sampled. We have developed observing and re¬ 
duction techniques for the GBT that result in high qual¬ 
ity spectral line data cubes (Boothroyd et al. 2011) with 
the requisite sensitivity for studies at low column density. 

An H I spectrum, whether for a single line of sight 
or averaged over a region, generally has emission spread 

^ The National Radio Astronomy Observatory is a facility of 
the National Science Foundation operated under cooperative agree¬ 
ment by Associated Universities, Inc. 
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over a range of frequency, which through the Doppler ef¬ 
fect is interpreted as radial velocity, in this paper tlsr 
relative to the Local Standard of Rest (LSR) and here¬ 
after denoted simply v. 

1.1. Insights from HI data 

This kinematic information in H I spectra provides es¬ 
sential diagnostics of various physical properties of the 
gas. H I gas in the local neighborhood is identifiable 
by its low |t|, while gas in the halo appears at both in¬ 
termediate and high velocities, for example in NEP in 
the ranges —80 < v < —20 kms“^ and —140 < v < 
—80 kms“^, respectively (for the GHIGLS fields, it turns 
out that the non-local gas is only at negative velocities.) 
Wakker & Woerden (2013) denote the latter two both as 
“high-velocity clouds.” However, for the three “compo¬ 
nents” that can be distinguished via velocity we prefer 
the standard terminology LVG, IVG, and HVC, respec¬ 
tively, because there is an underlying physical distinction 
between these components of interstellar gas, irrespec¬ 
tive of the apparent kinematics. There is a component 
of Galactic gas with an interesting history, circulating in 
a “Galactic fountain” (Shapiro & Field 1976; Bregman 
1980). Its distinctive motion projects into a radial veloc¬ 
ity that often sets it apart from local LVG emission, in 
the IVG range. On the other hand, there is gas that ap¬ 
pears to be of extragalactic origin (perhaps it has been in 
other galaxies) and is now accreting on the Galaxy and 
destined for interaction with Galactic gas in the halo and 
disk. Its radial velocity is often in the HVG range. Gom- 
plex C is one such example, with a low metallicity and 
high deuterium-to-hydrogen ratio that point to a non- 
Galactic origin (Tripp et al. 2003; Sembach et al. 2004). 

Surveys like GHIGLS can provide new insight into 
these three components, LVG, IVG, and HVG. For ex¬ 
ample, correlated dust emission corroborates this dis¬ 
tinction. Analysis with the GHIGLS data (Planck Col¬ 
laboration XXIV 2011) shows that the IVG gas has an 
emission signature from embedded dust with a dust-to¬ 
gas ratio comparable to the LVG, whereas for HVC there 
is no detectable dust emission signature, consistent with 
a dust-to-gas ratio at least as low as implied by the low 
metallicity. We note that depending on the geometry, gas 
with physical properties similar to IVG or HVC might 
appear within the LVG velocity range and so not be dis¬ 
tinguishable by its velocity (Wakker & Woerden 2013); 
some other clue is then needed such as, in the case of 
HVC-like gas, low H I-correlated dust emission (Planck 
Collaboration Int. XVH 2014). 

There is a long history of using the power spectrum 
of an image to quantify the statistical properties of in¬ 
tensity fluctuations and structural information, including 
H I (Crovisier & Dickey 1983) and thermal emission by 
dust (Gautier et al. 1992). This structure is linked to 
a turbulent cascade of kinetic energy. For H I Lazarian 
& Pogosyan (2000) showed that the power spectrum of 
a velocity channel is a complex mixture of velocity and 
density fluctuations. For the typical steep power spec¬ 
tra in the ISM, the three-dimensional spectral index of 
density can be obtained from power spectrum analysis 
at high spatial frequencies by averaging enough velocity 
channels that the brightness temperature fluctuations are 
dominated by density fluctuations; otherwise the power 
spectrum is too shallow by up to one in the power-law 


exponent. As pointed out and illustrated by Dickey et al. 
(2001) using H I observations in the inner Galactic plane, 
different regions of the ISM could have different statisti¬ 
cal properties projected on the sky because of different 
fractions of gas in the cold and warm phases of the neu¬ 
tral medium (CNM and WNM, respectively), different 
optical depth effects, and geometry. In Figure 10 from 
their study of turbulent molecular clouds Hennebelle & 
Falgarone (2012) provide a comprehensive summary of 
the power law exponents for different tracers of the ISM 
at different scales. Even for the same tracer, there is a 
considerable range. Gompared to fields near the Galactic 
plane or in molecular clouds, the GHIGLS fields at in¬ 
termediate Galactic latitude are relatively simple lines of 
sight, and yet they offer the opportunity to explore both 
environmental differences and the structural properties 
of LVG, IVG, and HVC components separately. 

The H I line profile, which might consist of a number of 
peaks with different centroid velocities, can be segmented 
by Gaussian decomposition (Hand 2000; Verschuur 2004) 
(this can be applied to H I absorption spectra as well, 
Roy et al. 2013). Gaussian decomposition methods of¬ 
fer the opportunity to differentiate between components 
with different line widths (Hand & Kalberla 2007). In 
combination with absorption-line studies (Dickey et al. 
2003; Heiles & Troland 2003), components with broad 
or narrow line widths have been found to arise from the 
WNM and CNM, respectively (Dickey & Lockman 1990; 
Wolflre et al. 2003; Hennebelle & Falgarone 2012). For 
the CNM the line width is larger than simply the ther¬ 
mal line width because of turbulent broadening. Such an 
analysis of H I data from surveys like GHIGLS can illu¬ 
minate theoretical numerical modelling of the phases in 
the ISM (e.g.. Saury et al. 2014), and vice versa, leading 
to a better understanding of the dynamical formation of 
CNM gas in a thermally bistable medium and the Ailing 
factors of the phases. 

1.2. Dust and Gas 

Observations of dust emission integrate over all dust 
along the line of sight, regardless of velocity. The strong 
dust-gas correlation seen empirically in projection at 
high Galactic latitude (Boulanger & Perault 1988) is con¬ 
sistent with a strong correlation spatially in three dimen¬ 
sions. The dust-gas correlation has been used extensively 
to infer the physical properties of dust in various environ¬ 
ments (Deni & Burton 1992; Boulanger et al. 1996; Jones 
et al. 1995). Reach et al. (1994, 1998) have argued that 
excess infrared emission relative to H I (an excess “dust 
emissivity”) toward brighter cloud structures signals a 
phase transition to H 2 , untraced by CO at those col¬ 
umn and volume densities. Subsequent all-sky observa¬ 
tions (Planck Collaboration XIX 2011) indicate that such 
“CO-dark” gas is an important component widespread in 
the Galaxy. 

The kinematic content of the H I spectra is important 
for decoding information in dust emission maps. Mor¬ 
phological spatial detail in maps of H I varies as a func¬ 
tion of velocity and so any dust closely correlated with a 
velocity component of the gas leaves a related morpho¬ 
logical imprint in the dust emission map. This approach 
has been used to show that there is dust of signiflcant 
emissivity associated with IVG gas (Martin et al. 1994). 

The selection of flelds targeted by GHIGLS enables ex- 
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ploration of the different kinematics and spatial distribu¬ 
tions of H I gas and the dust evolution at diverse stages 
of Galactic evolution. GHIGLS data have been used in 
combination with Planck data on thermal dust emission 
to find the emissivity, opacity, and temperature of dust 
associated with both LVC and IVC gas (Planck Collabo¬ 
ration XXIV 2011). Pinheiro Goncalves (2013); Pinheiro 
Goncalves et al. (2013) extend the analysis to emissivi- 
ties characterizing non-equilibrium dust emission in the 
mid-infrared IRAS bands (see also Pinheiro Goncalves, 
D. et al. 2015, in preparation). 

Miville-Deschenes et al. (2005) have reported thermal 
dust emission associated with the HVC gas in the Spitzer 
First Look Survey (FLS) field. (The GET Spectral Pro¬ 
cessor (SP) H I data used are described in Lockman & 
Condon 2005; they have been reprocessed here.) How¬ 
ever, this is challenging because the dust emissivity of 
this lower metallicity extragalactic gas has apparently 
been depressed by a lower dust to gas ratio and also 
potentially by a lower dust temperature. Furthermore, 
chance correlation of the foreground LVC and IVC dust 
emission with the cosmic infrared background fluctua¬ 
tions or “anisotropies” (thus CIEA), which so far have 
not been separated from the dust maps, is a significant 
source of systematic uncertainty and no HVC-correlated 
dust emission has been found in the GHIGLS fields 
(Planck Collaboration XXIV 2011). 

From another perspective, dust emission in the Galaxy 
is a significant foreground contamination of the CIEA 
signal. However, by exploiting the tight relation be¬ 
tween Galactic dust and gas emission at low column den¬ 
sities, H I-correlated dust emission can be removed from 
the infrared and submillimeter maps, as done for COEE 
(Arendt et al. 1998). In combination with Planck and 
Spitzer data, GHIGLS data have been used for “clean¬ 
ing” in this way to facilitate analysis of the residual CIEA 
(Planck Collaboration XVHI 2011; Penin et al. 2012; 
Planck Collaboration XXX 2014). 

Dust emission also contaminates measurements of the 
cosmic microwave background, for example compromis¬ 
ing detection of a polarized E-mode signal from inflation 
(Planck Collaboration Int. XXX 2014; EICEP2/Keck 
Array and Planck Collaborations 2015). The Galactic 
magnetic field tends to be oriented parallel to the elon¬ 
gation of filamentary dust structures (Planck Collabora¬ 
tion Int. XXXH 2014; Planck Collaboration Int. XXXV 
2015). Eilamentary gas structures in the cold neutral 
medium might also be useful in tracing the orientation 
of the magnetic field, another way in which H I data 
might contribute to a more complete understanding of 
the ISM both phenomenologically and physically. Clark 
et al. (2014) have presented observations of slender, lin¬ 
ear H I features in the diffuse ISM at high Galactic lati¬ 
tude and found them to be oriented along the interstellar 
magnetic field as probed by starlight polarization. 

The plan of the paper is as follows. The selection of 
the GHIGLS fields surveyed is presented in Section 2. 
Appendix A discusses data obtained using the GET SP 
while Appendix E discusses reprocessing of archival data 
overlapping one of these fields obtained using the GET 
ACS. Section 3 describes spectral-line mapping using the 
GET and the data reduction pipeline developed to pro¬ 
duce a data cube. Various ways of visualizing the data 


in a cube are reviewed briefly in Section 4 to illustrate 
the GHIGLS data. Separation of emission from distinct 
components of gas at different velocities, LVC, IVC, and 
HVC, is discussed in Section 5. Section 6 discusses prop¬ 
erties of maps of the line integral Hhi (proportional to 
column density Vri) for LVC, IVC, and HVC compo¬ 
nents. The effect of 21-cm line opacity on Vri is ad¬ 
dressed in Appendix C. Appendix D evaluates the un¬ 
certainties in the data, particularly as applied to maps 
of the column density Vrr To explore possible issues in 
data calibration and reduction, GET data in a few of the 
targeted GHIGLS fields are compared to data from a new 
generation of wide-area H I surveys, in the north EEHIS 
(Appendix E) and in the south GASS (Appendix E). Sec¬ 
tion 7 examines angular power spectra of the Vri maps 
of the three components. Gaussian decomposition of line 
profiles is explored in Section 8. The structure of the 
cold neutral medium and its relationship to the Galactic 
magnetic field are investigated in Section 9. Section 10 
summarizes our conclusions. 

2. THE GHIGLS FIELDS 

Table 1 presents the fields that we surveyed with the 
GET SP (Appendix A) and the GET ACS (since 2005), 
in order of increasing Galactic longitude. We also repro¬ 
cessed archival data obtained with the ACS for 09A079 
(Appendix E). Eor each field the table lists the center co¬ 
ordinates, adopted name, size, subfield layout where rel¬ 
evant (Section 3), the number of repeated observations, 
the noise def as measured in emission-free channels (Ap¬ 
pendix D.l) per 0.8 kms“^ or 1.0 kms“^ channel (Sec¬ 
tion 3.1), and the scan orientation (Section 3.2), along 
lines of constant Galactic latitude (Galactic) or Declina¬ 
tion (ICRS). 

The raw H I spectral line data using the GET SP 
are available under the following proposal numbers: 
GET/02A-007, 02A-023, 02A-031, and 03E-030. Those 
using the GET ACS are under GBT/05C-009, 05C-021, 
06E-030, 06C-032, 07A-104, 08A-083, 08E-038, 09A- 
079, 09E-042, lOA-012, and lOA-078. The field names 
adopted here are recognizable designations for the data 
in the NRAO archive,^ except for DEN (both CDEN and 
HDEN), SPIDER (DDI), UMIN (NEP42), and 09A079 
for the archival data described in Appendix B. 

Eigure 1 shows the fields in the northern Galactic in¬ 
termediate latitude sky and inserts for three southern 
fields, including the one ACS field in the southern sky, 
MC, which covers part of the Magellanic Stream. The 
North Celestial Pole (NCP) is marked by an “x” and 
those fields scanned in Equatorial coordinates can be 
identified by their different orientation. A “-h” marks 
the North Ecliptic Pole (NEP). Although many fields 
are close neighbors, even intentionally overlapping, given 
their size a considerable variety of physical conditions has 
been probed, sometimes serendipitously but particularly 
because the fields were chosen with certain science goals 
in mind. 

2.1. Science Goals 

Below we describe some motivations for targeting par¬ 
ticular fields for our H I surveys, building on the gen- 
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Table 1 

The Constituent Fields in GHIGLS 


Field 

Name 

Size 

Subfield 

Layout 

Repeats 

[mK] 

Scan and 
Goordinates 

G056.98-81.50 

MG 

6° X 5° 


2 

81 

Galactic 

G058.10+68.55 

EOOTES 

12? 4 X 4° 

5x1 

5(2) 2,3 

3/ 1(2) 5 

70(110) 4 

50/ 80(60) 4 

IGRS 

G067.74+67.73 

Necklace 

2° X 2° 


2 

73 

Galactic 

G071.00+41.47 

0X3 

1?9 X 2?4 


1 2 

90 4 

Galactic 

G085.33+44.28 

N1 8 

5° X 5° 


2 7 

71 

Galactic 

G087.95+59.05 

G86 

5° X 5° 


3 

59 

IGRS 

G088.32+34.89 

FLS 

3° X 3° 


82 

82 

IGRS 

G091.38+47.95 

MRK290 

4° X 4° 


1 

95 

Galactic 

G092.24+38.43 

DRAGO 

5° X 5° 


3 

60 

IGRS 

G096.27+59.91 

GROTH 

2° X 2° 


32 

77 4 

IGRS 

G096.40+30.03 

NEP 

12° X 12° 

3x3 

3 

64 

Galactic 

G115.62+30.40 

UMIN 

4° X 3° 


1 

100 

Galactic 

G125.00+27.42 

POL 

6° X 10° 


1 

111 

Galactic 

G125.03+37.36 

POLNOR 

6° X 10° 


2(1) 8 

75(110) 

Galactic 

G125.37+41.67 

MRK205 

4° X 4° 


1 

104 

Galactic 

G125.89+54.84 

DFN 

2° X 2° 


5 2 

45 ^ 

IGRS 

G132.37+47.50 

SP 

5° X 5° 


2 

71 

Galactic 

G134.95+54.13 

UMl 

4° X 4° 


1 2 

83 4 

IGRS 

G134.98+39.97 

SPIDER 6 

10° X 10° 

5x5 

2/ 1 9 

75/ 105 

Galactic 

G135.36+30.29 

SPG 

12° X 9?5 

3x2 

1/ 2 10 

103/ 72 

Galactic 

G143.94+28.02 

1H0717 

4° X 4° 


1 

97 

Galactic 

G144.25+38.56 

UMA 

9° X 9° 

3x3 

1 

107 

Galactic 

G145.68+23.35 

HS0624 

4° X 4° 


1 

101 

Galactic 

G147.46+44.09 

UM3 

4° X 4° 


1 2 

82 4 

IGRS 

G148.65+52.21 

09A079 

5° X 4° 

irregular 

1 to 4 

62 

Galactic 

G152.31+53.31 

UM2M 

4° X 4° 


1 2 

86 4 

IGRS 

G152.44+25.63 

MS0700 

4° X 4° 


1 

101 

Galactic 

G155.76+37.00 

UMAEAST 

10? 5 X 6° 

3x1 

1 

107 

Galactic 

G156.38+32.57 

LOOP4 

4° X 4° 


1 

122 

Galactic 

G156.45+54.06 

NGG3310 

4° X 4° 


1 

95 

Galactic 

G158.32+28.75 

MRK9 

4° X 4° 


1 

96 

Galactic 

G164.84+65.50 

AG 

5° X 5° 


2 

71 

Galactic 

G170.02-59.90 

SUEA 

2° X 2° 


22 

83 4 

IGRS 

G172.01+26.84 

MEM23 

4° X 4° 


1 

99 

Galactic 

G175.36+43.38 

091346A 

4° X 2° 


1 2 

83 4 

IGRS 

G179.50+65.03 

MRK421 

4° X 4° 


1 

95 

Galactic 

G223.57-54.44 

GDFS 

2° X 2° 


52 

93 4 

IGRS 


^ Average noise in emission-free channels of width 0.8 kms“^ (except 1.0 kms“^ where noted"^). 

^ Using the GET SP (Appendix A). 

^ Central 4° x 4° subfield with two repeats. Inner 3?5 x 3?5 three further repeats, totalling five. 

Channel spacing 1.0 kms“^. 

^ Four overlapping 2?5 x 4° subfields fianking the central GET SP subfield symmetrically. Three repeats, 
except western-most only once with its central 0.24 fraction in S repeated twice. All regridded to 1.0 kms“^ 
channel spacing of GET SP data. 

® A DRAG Deep Field (Elagrave, K. et al. 2015, in preparation). 

^ From data taken at times that minimize the stray radiation correction (Appendix D.5). 

® Upper 0.53 fraction of field repeated twice. 

^ Two repeats on the inner 3 X 3 of the 5x5 subfields. 

Three 4° x 6° subfields and above these three 4° x 2?5 subfields; only the western pair of subfields 
observed twice. Additional 1° strips on north and west to provide an overlap with SPIDER and with POL 
and POLNOR, respectively. 

Reprocessed archival GET ACS data toward the Lockman Hole (Appendix E). 


eral discussion in Section 1. Abstracts of our propos¬ 
als making use of the GBT ACS are also available at 
library.nrao.edu/proposals. 

However, it can be noted, certainly in retrospect, that 
data obtained in a focused proposal are often useful 
for addressing several of the distinct science goals high¬ 
lighted in others. This is illustrated by the data for 
many GHIGLS fields that proved useful for the above- 
mentioned Planck studies of the GIB A and Galactic dust. 

2.1.1. Extragalactic Windows 

H I surveys are used first to identify regions of low 
Galactic column density and then to characterize or re¬ 


move foregrounds that otherwise compromise extragalac¬ 
tic science. With the GBT SP, the foreground H I has 
been explored in a number of extragalactic or “cosmic” 
windows for deep multiwavelength observations of the 
distant universe. The fields UM2M together with 09A079 
and NGC3310 (both with the GBT ACS) overlap the 
Lockman Hole, the region with the lowest integrated H I 
emission (Figure 1). The Lockman Hole has been studied 
previously in H I at 10' resolution using the NRAO 300-ft 
telescope (Lockman et al. 1986; Jahoda et al. 1990) and 
was an important field for the Spitzer SWIRE survey 
(Lonsdale et al. 2003) and many other studies, e.g., with 
XMM-Newton (Hasinger et al. 2001). UMl has similarly 
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Figure 1. Locations of GHIGLS fields. Fields are rectangular in the local GLS projection appropriate to the GTE mapping but in 
this stereographic (STG) projection the outlines are slightly distorted. Background image is the velocity-integrated H I emission from the 
all-sky LAB survey (Kalberla et al. 2005), expressed as column density Nm in units 10^® cm“^. The North Gelestial Pole (NGP) is marked 
by an “x” at 122?93, 27? 13. The GBT fields cover the entire NGP Loop and its interior. The North Ecliptic Pole (NEP) is marked by a 
at 96?38, 29?81. MG and two other southern fields are shown to scale in the insets on GLS grids. 


low LVC emission, but has much brighter IVC emission 
and some HVC as well. 

In addition to those discussed further below, BOOTES 
(extended with the ACS) and N1 (with the ACS), 
GHIGLS GBT SP fields span many other notable 
windows: 0X3 - the Hercules field in the Oxford- 
Dartmouth Thirty Degree Survey (MacDonald et al. 
2004); GROTH - the All-Wavelength Extended Groth 
Strip International Survey (AEGIS) region (Davis et al. 
2007); DEN - Chandra Deep Eield North/Hubble Deep 
Eield North/GOODS-N (Hornschemeier et al. 2001); 
SUBA - the Subaru/XMM-Aei(;^6>n Deep Survey Eield 
(Ouchi et al. 2001) (also SWIRE); and CDES - Chandra 
Deep Eield South/GOODS-S (Giacconi et al. 2001) (also 
SWIRE). We note that these areas include four of the 
five (tiny) CANDELS fields (Grogin et al. 2011). 

2.1.2. High Velocity Infall 

The European Large-Area ISO Survey (ELAIS) N1 
field is an extragalactic window also targeted by SWIRE 
whose low column density can be appreciated in Eigure 1. 
The complex H I structure has been studied at 1' resolu¬ 
tion using the Synthesis Telescope (ST) at the Dominion 
Radio Astrophysical Observatory (DRAO^), mosaicing 
data from 76 closely-packed pointings to produce a sensi¬ 
tive “DRAO Deep Eield” for which our GBT ACS data, 
the GHIGLS N1 field, provide the short spacing infor¬ 
mation (Blagrave, K. et al. 2015, in preparation). Unlike 
the Lockman Hole, this region has striking HVC emis- 

® WWW.nrc-cnrc.gc.ca/eng/solutions/facilities/drao.html 


sion, enabling a search for HVC-correlated dust emission. 
The HVC gas morphology is not immediately obviously 
traced by dust, and in fact it turns out quantitatively, at 
least at the GBT angular resolution, that the HVC has a 
low dust emissivity (Planck Collaboration XXIV 2011). 

Along with N1 and ELS, we extended our study of 
dust emissivity in HVC complexes by selecting regions 
in which the HVC to IVC plus LVC column density 
in the LAB survey showed a high contrast: MC, SP, 
and AG. (There are also significant knots of high con¬ 
trast HVC emission contained within SPC, UMA, and 
UMAEAST.) Their low total column density makes these 
fields suitable for CIBA studies as well (Planck Collab¬ 
oration XVHI 2011; Penin et al. 2012; Planck Collab¬ 
oration XXX 2014), but we found that the CIBA sig¬ 
nificantly contaminates studies of (HVC) H I-correlated 
dust emission (Planck Collaboration XXIV 2011). 

2.1.3. The Galactic Fountain 

DRACO (Herbstmeier et al. 1993) and G86 (Martin 
et al. 1994) were selected because of their prominent 
and distinctive IVC gas which has a clear dust signa¬ 
ture. Targeted Herschel dust emission and DRAG ST 
H I observations were carried out for these two fields as 
well. Such fields are ideal for searching for differences in 
dust properties between local and IVC gas, evidence of 
dust evolution in different IVC environments. The dy¬ 
namics is fascinating and in DRACO the transition from 
atomic to molecular gas can be studied. 

The central region of our BOOTES field targeted the 
extragalactic window known as the NO AO Deep Wide 
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Field Survey (Jannuzi & Dey 1999) and the focus of a 
Spitzer MIPS wide/shallow survey (Dole et al. 2004). It 
is actually crossed by a band of Galactic dust emission 
within which we subsequently discovered anomalously 
high dust emissivity in one faint component of IVC gas 
(Lockman et al. 2005); to explore this further, with the 
GET ACS, we added flanking subfields along the direc¬ 
tion of the band. In IRAS 60 and 100 /rm dust maps 
of the Bootes region there appears a remarkable feature 
that we dubbed the Necklace and for which we carried 
out a small H I survey to study the dust emissivity. 

2.1.4. North Ecliptic Pole Foreground 

We mapped H I in a large region centered on the NEP. 
As far as gas and dust content is concerned, there is 
nothing special a priori about this intermediate latitude 
location. However, for many all-sky surveys by satellites 
that scan along great circles passing near the ecliptic 
poles, the NEP region is special in having more coverage 
and hence data products with lower uncertainty. This is 
the case for Einstein^ IRAS, COBE, ROSAT, AKARI, 
and recently Planck^ as illustrated in the coverage map 
(see, e.g., figure 5 in Planck Collaboration I 2014) and 
related variance maps. Thus the NEP region is the focus 
for many studies, e.g., recently with Chandra (Krumpe 
et al. 2015). Eor several science applications using these 
data it is important to understand the foreground from 
the Galaxy. 

H I in this region was surveyed previously using the 
140-ft telescope at NRAO (Lockman et al. 1986). Our 
NEP GBT ACS data cover a larger region with higher 
angular resolution and much reduced noise and system¬ 
atic uncertainty. The observations, repeated three times, 
have a good signal to noise ratio, there is significant LVC, 
IVC, and even HVC emission, and the field is large, all 
reasons why we chose NEP data for illustrations in this 
paper. 

2.1.5. The North Celestial Pole Loop 

The first large field targeted with the GBT ACS was 
SPIDER, at the top of the arch of the North Celestial 
Pole Loop (NCPL), a giant gas structure north of the 
Galactic Plane with a cylindrical morphology (see Mey- 
erdierks et al. 1991). In SPIDER the complex LVC H I 
structure has been studied with the DRAG ST, mosaic- 
ing data from 91 pointings to produce a second DRAG 
Deep Eield. We have also mapped the stunning dust 
emission there at even higher resolution with Herschel. 

Ultimately we explored H I throughout the spider’s 
web, along the entire NCPL extending from PGL to 
UMAEAST, and also H I in its “interior” below the arch 
with SPC and several fields toward higher longitude. See 
Eigure 1. PGL covers the Polaris flare, which is actually 
dominated by molecular gas, not H 1 . Herschel dust emis¬ 
sion maps and DRAG ST H I observations were acquired 
for parts of PGL and UMA (Blagrave, K. et al. 2015, in 
preparation). 

2.1.6. Molecular Hydrogen 

In some of the GHIGLS fields analyzed in Planck Col¬ 
laboration XXIV (2011) there was evidence from excess 
submillimeter dust emission for gas not traced by H I, 
suggesting the presence of H 2 . To explore evolutionary 


aspects of this phenomenon in relation to different envi¬ 
ronmental and kinematic factors at intermediate Galac¬ 
tic latitude, we focused on fields containing an AGN for 
which EUSE observations have been used to character¬ 
ize H 2 (Gillmon et al. 2006; Wakker 2006) and one field 
containing a small, non star forming molecular cloud, 
MBM 23 (Magnani et al. 1985). Many of these fields 
have been observed with the DRAG ST as well, for which 
our GBT observations provide the short spacing data. 

3. SPECTRAL-LINE MAPPING WITH THE GBT 

Surveys of H I spectra in sizeable fields were carried out 
using on-the-fly (GTE) mapping. As described in detail 
by Boothroyd et al. (2011), data were acquired with the 
GBT ACS by in-band frequency switching, with spec¬ 
tra recorded every 4 s in two independent polarizations. 
The velocity coverage was —1700 < v < +900 kms“^ 
in independent 0.161 kms“^ channels, from which was 
extracted the central —450 <v< 400 kms“^. With the 
broad spectral coverage, emission from LVC, IVC, and 
HVC gas is all accessible. 

In our H I surveys, scans of the sky up to 6° long were 
made. Eor practical purposes larger fields were broken 
up into smaller subfields with scan dimension between 
2° and 4° that were mapped separately. The scans are 
still sufficiently long that the choice of GTE mapping 
greatly reduces the system overhead. The mapping speed 
is typically one deg^ per hour. As described below, we 
mapped over 1000 deg^ using the GBT ACS. Gbserva- 
tions of many fields were repeated both to increase the 
sensitivity and to examine the reproducibility of the data, 
leading to better maps. 

Aspects of the mapping, processing of the spectra, and 
data cubes are discussed in this section. The accuracy of 
the data has been discussed by Boothroyd et al. (2011) 
and this is complemented by the analysis in Appendix D. 

3.1. Processing of the Spectral Data 

The desired quantity to be measured, designated sim¬ 
ply by Tb in this paper, is the H I brightness temper¬ 
ature averaged over the main beam. As described by 
Boothroyd et al. (2011), the data reduction process in¬ 
volves removal of radio frequency interference (REI), pro¬ 
ducing the frequency-switched spectra including a cubic 
spline interpolation of the data onto a common grid in 
u, spectral smoothing and resampling to 0.807 kms“^ 
channels^^ (or 1.03 kms“^), averaging the two polariza¬ 
tions, calibration of the intensity scale to antenna tem¬ 
perature, Ta, calculation and subtraction of the stray ra¬ 
diation spectrum, correction for the main beam efficiency 
and the atmosphere, interpolation of the sampled spec¬ 
tra from all repeated observations into a data cube, and 
finally removal of a small (low-order and low-amplitude) 
instrumental baseline for each interpolated spectrum in 
the data cube. In some cases, the presence of a galaxy 
makes it impossible to fit a reliable baseline and thus 
no baseline is removed. A mask developed for the data 
product release reflects these cases. 

The repeated observations were also used to pro¬ 
duce independent data cubes for diagnostic analysis (see 

We filtered the spectra with an eleven-channel Hanning kernel 
(nine non-zero weights) and then sampled every fifth channel to 
produce independent data. The effective resolution is 1.0 kms“^. 
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Boothroyd et al. 2011 and Appendix D). Similarly, we 
have compared data products obtained by processing 
the two independent polarizations separately (Boothroyd 
et al. 2011) to check for reliability. 

3.2. Mapping, Gridding, and Resulting Angular 
Resolution 

The native primary beam of the GBT is only mildly el¬ 
liptical with full width half maximum (FWHM) 9.T and 
9.'0 in the cross-elevation (azimuthal) and elevation direc¬ 
tion, respectively. The integration time (4 s) and tele¬ 
scope scan rate were chosen to sample the spectrum every 
3.'5 in the in-scan direction, more finely than the Nyquist 
interval, 3.'86. Beam broadening in OTF mapping can be 
made negligible by sampling at a smaller fraction of the 
Nyquist rate (Mangum et al. 2007), ideally at half that 
interval. Our actual sampling amounts to a convolution 
of the beam with a 3.'5 boxcar, extending the beam to 
an effective FWHM of 9.'4 in the in-scan direction (not 
aligned along either azimuth or elevation). 

Fields were mapped boustrophedonically, scanning the 
telescope in one direction (Galactic longitude or Right 
Ascension), with steps of 3.'5 in the orthogonal (cross¬ 
scan) coordinate direction before the subsequent reverse 
scan. The quality of the spectra was monitored closely 
during most sessions. Occasionally artifacts from the 
spectrometer would compromise a spectrum and it was 
flagged; usually it was possible to schedule an observa¬ 
tion to replace the entire scan. 

The individual spectra were gridded into data cubes. 
For these large fields, the equal-area Global Sinusoidal 
(GLS) coordinate projection (Calabretta & Greisen 
2002) was chosen to complement our equal-area sam¬ 
pling strategy. GLS is now replaced by the Sanson- 
Flamsteed (SFL) projection, which can be considered 
equivalent as used here. Gridding data from a large 
region to another projection, like TAN, would produce 
pixels (spectra in the cube) with non-uniform coverage. 
The data are mapped to the grid following a convolu¬ 
tion. The simplest convolution would be a pillbox, the 
naive-mapmaking strategy of assigning each sample to 
its nearest neighboring grid cell. With GLS projection, 
a grid can be defined with 3.'5 pixels corresponding to 
where the data were sampled, in which case each pixel 
contains the actual spectrum recorded. This is very use¬ 
ful for diagnostic purposes as we assess the accumulation 
of error in individual spectra. 

However, pillbox is not the ideal gridding for our data. 
Instead we used an optimal tapered Bessel function for 
interpolation (Mangum et al. 2007), available in Classic 
AIPS with the task SDGRD. This approach also allowed 
us to account for small telescope pointing deviations in 
the in-scan direction relative to the expected uniform 
raster. Convolution with a Bessel function, being the in¬ 
verse Fourier transform of a top hat, is the equivalent of 
equally sampling all spatial scales in the Fourier domain 
with a two-dimensional boxcar function of size extending 
to the resolution limit D/X = 1/7.25 arcmin“^. It pre¬ 
serves the true power on small spatial scales while avoid¬ 
ing the introduction of noise on scales smaller than the 
beam. The corresponding “width” of the Bessel function 
is 3.'75. Because the observed fields are finite (i.e., no in¬ 
formation exists beyond the edge of the observed field), 
the Bessel function is tapered by a Gaussian. The sup¬ 


port size (the radius beyond which the tapered Bessel 
function is truncated) was chosen to correspond to a 
zero intercept of the Bessel function. Both the first 
and second of these intercepts, corresponding to 7.'5 and 
11.'25, were considered. A larger support size allows for a 
broader Gaussian taper resulting in less modification of 
the beam profile; while this results in loss of information 
at the map edge, for our large maps this loss is inconse¬ 
quential. Thus a support size of 11.'25 was selected along 
with a Gaussian taper “width” (as defined in the AIPS 
task SDGRD) of 9.'16, which corresponds to a Gaussian 
with FWHM = 15.'26. The width of the Gaussian taper 
is chosen to be as broad as possible, while also main¬ 
taining a smooth transition from the (otherwise infinite 
in nature) Bessel function to a value of zero beyond its 
support size (i.e., its truncation radius). 

As is the case with any non-pillbox grid system, the 
beam ends up being broadened slightly in both in-scan 
and cross-scan directions. The resulting effective beam 
for our modified Bessel function gridded data cubes can 
be approximated by an elliptical Gaussian of size 9.'55 x 
9.'24 (FWHM). 

Data cubes were constructed in Classic AIPS, combin¬ 
ing spectra from the two polarizations and from repeated 
observations as appropriate. A weight map was created 
using SDGRD during convolution of the spectra into the 
GLS grid and this map was used to remove a few pixels 
along the edge of the map, leaving a cube with uniform 
coverage. Spectra from the subfields mapped were com¬ 
bined directly into a common GLS grid (with the same 
3.'5 pixel size), providing a mosaic of the entire region. 
Because the areas covered by the subfields are rotated 
slightly compared to this common grid, there is a saw¬ 
tooth pattern to the weighting function along the mosaic 
edges. This too was removed for the final data prod¬ 
ucts. These cropped mosaics have near-uniform coverage 
(weight), increasing where the subfields overlap. 

These cubes were examined for any evidence of the ef¬ 
fects of an anomalous spectrum that might have been 
missed in the previous flagging, in which case the bad 
spectrum was identified and flagged and the cube re¬ 
made, where possible including a replacement scan; any 
remaining locations affected by missing/flagged spectra 
are identifiable in the weight and noise maps. Finally 
a baseline was fit to the gridded data in each pixel and 
subtracted. A mask was developed to record the few pix¬ 
els for which a baseline fit to the spectrum was deemed 
not feasible, usually due to the presence of emission from 
a galaxy; the uncorrected spectra are still very good be¬ 
cause of the intrinsically flat baselines of the GBT ACS. 

On the archive of GHIGLS data,^^ for each field there 
is a FITS file with extensions as follows: 0, the cube 
of spectra Tb in K; 1, a mask recording the few pixels 
with no baseline fit removed (many fields do not require 
this extension); 2, a map of the noise def as measured in 
emission-free channels of each spectrum, in K, and 3, the 
(relative) weight map. In the noise map can be seen the 
reduction in noise from overlapping/repeat observations, 
the subfield layout where relevant (e.g.. Figure 15), and 
the rare increases in noise because of flagged spectra. 
These are a direct consequence of the weighting and so 
appear in the weight map as well. 

www.cita.utoronto.ca/GHIGLS 
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3.3. Mosaiced Data 

The data for the tabulated GHIGLS fields have been 
analysed separately, but to explore larger-scale connec¬ 
tions between features in individual fields we have com¬ 
bined the data in some adjacent/overlapping fields into 
mosaics, in the same way that subfields were combined 
for a given target. In this reprocessing, we first ensured 
that the individual spectra had a common velocity grid. 
Then we executed the last two steps of our pipeline, in¬ 
terpolation of the spectra onto a common GLS grid, and 
removal of a small instrumental baseline fit at each pixel. 
FITS files of the mosaics are available on the GHIGLS 
archive. 

We have made a small mosaic LHM covering the low 
column density area in the Lockman Hole region by com¬ 
bining the data from NGG3310, UM2M, and 09A079. 
Because of the overlap the sensitivity is improved over 
the central region, as can be seen in the noise map. 
Glose inspection shows that of the three fields, UM2M 
has slight residuals from the stray radiation correction. 

We have combined data for NEP, ELS, and DRAGO 
into a mosaic named NED. This reprocessing also in¬ 
cludes data for a small 3 deg^ map KnotN under proposal 
GBT/06B-030 exploring the IVG and an embedded 16 
deg^ field H1821 under GBT/09G-042, not documented 
in Table 1 or Eigure 1 but evidenced in the noise map. In 
this mosaic it can be seen how the distinctive LVG fila¬ 
ment in ELS extends into DRAGO and is roughly parallel 
to the filaments already mentioned in NEP. 

As Eigure 1 shows, GHIGLS has extended cover¬ 
age of the NGPL and so we have made a large mo¬ 
saic thus named. This reprocessing also includes data 
from a few embedded surveys, namely Knot A, B, and 
G under GBT/05G-021, our exploratory focus on high- 
contrast HVG emission, plus fields called PG0804 under 
GBT/09B-042 and LISZTA, LOOPl+2, and LOOP4B 
under GBT/lOA-012. This shows the spectacularly com¬ 
plex structures in the LVG emission collectively defining 
the arch, but also fascinating IVG and HVG emission as 
well. 

4. VISUALIZATION OF A DATA CUBE 

The complex H I spectral information Tb(x,^,u) as¬ 
sembled in a data cube can be visualized in a number 
of ways. The renderingin Eigure 2 for the NEP field 
reveals the characteristic clumpy structure of Tb in the 
three dimensions of the cube. We note, however, that 
V does not necessarily map directly into a third spatial 
dimension. The sense of connectedness and separation of 
the structures can be reinforced by interactively chang¬ 
ing the viewing angle (azimuth and elevation). We have 
created a movie to illustrate this, smoothly varying the 
viewing angle so that the cube appears to tumble before 
eventually returning to the frame shown in Eigure 2. 

Eigure 3, also for NEP, shows three complementary 
slices of the cube: “a channel map” and two “position- 
velocity diagrams.” Position-velocity diagrams reveal not 
only the complex structure in the gas, but also regularity 
such as distinct ranges in v in which the emission is con¬ 
centrated (e.g., distinct LVG, IVG, and HVG are clearly 

Created with SAOImage DS9 (http://ds9.si.edu) using a 
simple ray-trace algorithm with the Maximum Intensity Projection 
method. 



Figure 2. Rendering of the Tb(x,y, n) data in the spec¬ 
tral cube for the NEP field between —210 and 60 kms“^. 

The colorbar is linear in Tb (K). Changing the viewing 
angle can reinforce the sense of connectedness and sep¬ 
aration of the structures, as illustrated in a movie of a 
tumbling cube beginning and ending at the frame shown here 
(www.cita.utoronto.ca/GHIGLS/MOVIES/TUMBLE/GHIGLS_NEP_tumble.mp4^ ). 

^All movies are also available in “ogg” format. 

seen in the {i, v) diagram shown at the lower right). Gon- 
nectivity in these diagrams, or more generally in the ren¬ 
dered cube, can suggest a physical relationship between 
the different velocity components. A sequence of chan¬ 
nel maps (or orthogonal planes) can be combined as the 
frames of a movie file. Eor the NEP we have made three 
such movies, within which the images in Eigure 3 are sin¬ 
gle frames (see caption). Eor the other GHIGLS fields, 
equivalent movies including tumbling cubes can be found 
on the GHIGLS archive. 

The individual spectra at the lower left in Eigure 3 can 
be compared to the mean spectrum for the cube in Eig¬ 
ure 4; across the field the spectrum changes dramatically 
in the LVG, IVG, and HVG ranges. 

A striking phenomenon seen over the LVG range of the 
NEP cube is a series of filamentary structures all running 
roughly diagonally across the field. Some of these can be 
seen in the single channel map in Eigure 3. As discussed 
in Section 9, these appear to be aligned with the local 
Galactic magnetic field. 

5. SEPARATING LOW, INTERMEDIATE, AND HIGH 
VELOGITY GAS 

Separating the H I emission into contributions from 
LVG, IVG, and HVG gas is challenging when the velocity 
ranges of the components overlap. Eor NEP the HVG 
emission is relatively weak but in Eigures 2 and 3, and 
even in the average spectrum in Eigure 4, it is fairly well 
separated. The IVG emission is not as strong as LVG 
and while the existence of IVG is apparent in the average 
spectrum in Eigure 4 it is not obvious exactly where to 
make the separation from LVG. 

A useful diagnostic is based on the standard deviation 
about the mean of a channel map, because it depends not 
just on the presence of signal but also on the fluctuations. 

Thus the standard deviation spectrum takes advantage of 
the rich structure within the cube and so when the emis¬ 
sion corresponding to LVG, IVG, and HVG components 
is not immediately distinguishable in the mean spectrum 
it offers the potential to separate the cube into distinct 
velocity ranges (Planck Gollaboration XXIV 2011). 
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Figure 3. Maps of on three complementary orthogonal 
slices of the H I spectral cube for the NEP field. Upper 
right: Channel map (^,6) ai v = —4.0 kms“^, viewed 

from the high v face of the cube as rendered in Figure 2. 
This is a single frame from a movie running through chan¬ 
nel maps of the entire cube, starting from the highest v 
(www.cita.utoronto.ca/GHIGLS/MOVIES/LB/GHIGLS_NEP.mp4) . 
Striking filamentary structures cross the field diagonally, roughly 
parallel to the direction of the Galactic magnetic field (Section 9). 
Solid vertical or horizontal lines represent the cross-section 
cut from which the adjacent complementary position-velocity 
diagrams are built. The two position-velocity diagrams are 
as follows. Lower right: {i,v) at fixed b = 25?4; the view is 
from the low-6 face of the cube in Figure 2 with the v axis 
foreshortened. This is a single frame from a movie running 
through (^, v) slices of the entire cube, starting from the lowest 6 
(www.cita.utoronto.ca/GHIGLS/MOVIES/LV/GHIGLS_NEP_lv.mp4) . 
Vertical dashed lines mark the longitudes along which profiles 
of the intensity in this diagram are shown in the panel at the 
lower left; these profiles are the spectra at two of the pixels 
marked by bulls-eyes in the channel map. Upper left: (v, 6) for 
i ~ 94°, more precisely the plane through {i,b) = (94? 1,25? 4); 
the view is from the high-^ face of the cube in Figure 2 with the 
V axis foreshortened. This is a single frame from a movie running 
through (i;, 6) slices of the entire cube, starting from the highest £ 
(www.cita.utoronto.ca/GHIGLS/MOVIES/VB/GHIGLS_NEP_vb.mp4) . 
Horizontal dashed lines mark the latitudes of the spectra shown. 
Lower left: Distinctive individual spectra for the three pixels 
marked by the bulls-eyes in the channel map, along the dashed 
lines in the position-velocity diagrams. The spectrum at the 
special pixel (94? 1,25? 4) can be seen in both position-velocity 
diagrams. 


We have evaluated this approach using two cubes 
from the hydrodynamical simulations of the structure 
of the thermally bistable and mildy turbulent atomic 
gas in the local ISM (Saury et al. 2014). These in¬ 
clude turbulent forcing with a mixture of compressive 
and solenoidal modes whose partition is set by the spec¬ 
tral weight ( (ranging from 0 for compressive to 1 for 
solenoidal), a turbulent forcing amplitude and ini¬ 
tial density no, and for investigating the approach to 
component separation here, the exact choice of param¬ 
eters/simulations does not matter. Each of these cubes 
consists only of LVC. However, we combined them by 
summing Tb(x, y, v) after translating one cube by Aughift, 
thus producing an IVC component as well. We added 



-200 -150 -100 -50 0 50 
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Figure 4. Spectral profiles of the mean (solid) and median 
(dashed) values of Tb for the NEP field (left axis) and of the stan¬ 
dard deviation about the mean of Tb('^) (dotted, ctt, right axis). 
The complexity of H I structure along the line of sight compli¬ 
cates the separation of the different radial velocity components. 
Although more blended in the mean and median profiles, the LVC, 
IVC, and HVC components become more distinct in the standard 
deviation spectrum. 

noise characteristic of observations of the GET fields. 
When Aughift exceeds the typical FWHM of the line 
profile, the separation of components is unambiguous, 
whether using the mean or standard deviation spectrum. 
But the separation becomes more challenging for smaller 
Aughift and the superior utility of the standard deviation 
spectrum is clear. We quantified this by comparing the 
column densities (Section 6) computed for the two ranges 
to the actual column densities in the original cubes. For 
conditions like encountered in the actual GHIGLS fields, 
the internal boundaries between LVC and IVC and be¬ 
tween IVC and HVC can be determined to better than 
a few kms“^. 

The application to NEP is shown in Figure 4, from 
which the selected bounding velocities are [47.9, —20.5, 
-83.3,-199.9] kms“^. Table 2 shows the adopted ve¬ 
locities for component boundaries in the GHIGLS fields. 
These values correspond precisely to particular channels 
boundaries in the data cube, but depending on the field 
the actual internal boundaries can be uncertain by a few 
channels. 

If the separation is ambiguous, it becomes a potential 
systematic error in whatever dependent analysis is be¬ 
ing carried out and so a sensitivity analysis needs to be 
done. Consider the case of distinct dust emissivities as¬ 
sociated with LVC and IVC components. Our approach 
has been to correlate simultaneously the LVC and IVC 
column density with a dust emission map from the all-sky 
Planck survey (Planck Collaboration XI 2014), with the 
two emissivities as free parameters. When the compo¬ 
nent separation seemed somewhat arbitrary, the velocity 
boundary between the LVC and IVC components was 
changed incrementally over the range of possible veloci¬ 
ties, producing a pair of emissivities and an rms of the 
residual dust map for each case. The velocity bound¬ 
ary corresponding to the minimum value of the rms was 
consistent with that determined independently from the 
standard deviation spectrum. The derived values of the 
emissivities were robust and the uncertainty of the ve¬ 
locity boundary did not significantly increase the uncer- 
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Table 2 

Component Velocity Boundaries (km s“^) 


Name 

LVG 

IVG 

HVG 

MC 

50.4 

-6.0 

-39.8 

-240.1 

BOOTES 

29.5 

-17.9 

-69.4 

-111.6 

Necklace 

14.9 

-9.2 

-50.3 


0X3 

50.2 

-28.2 

-95.1 

-179.6 

N1 

57.6 

-10.8 

-59.9 

-151.6 

G86 

23.8 

-26.9 

-64.8 

-137.2 

ELS 

27.8 

-19.7 

-114.6 

-210.3 

MRK290 

35.1 

-47.9 

-80.9 

-175.0 

DRAGO 

23.8 

-8.4 

-72.8 

-209.5 

GROTH 

27.5 

-37.4 

-66.3 

-125.0 

NEP 

47.9 

-20.5 

-83.3 

-199.9 

UMIN 

50.4 

-22.9 

-58.3 

-129.9 

POL 

31.8 

-48.7 

-97.0 

-140.4 

POLNOR 

31.8 

-26.9 

-97.0 

-220.0 

MRK205 

39.9 

-35.0 

-100.2 

-220.0 

DEN 

50.2 

-28.2 

-105.4 

-165.2 

SP 

31.8 

-30.2 

-105.0 

-190.2 

UMl 

60.5 

-24.1 

-93.1 

-175.4 

SPIDER 

39.9 

-14.9 

-88.1 

-159.7 

SPG 

32.6 

-47.9 

-100.2 

-216.8 

1H0717 

50.4 

-36.6 

-109.8 

-210.3 

UMA 

50.4 

-26.1 

-84.9 

-193.5 

HS0624 

39.9 

-33.4 

-68.0 

-150.0 

UM3 

29.5 

-25.1 

-84.8 

-170.3 

09A079 

60.0 

-30.2 

-80.1 

-159.7 

UM2M 

39.8 

-30.2 

-79.7 

-125.0 

MS0700 

36.7 

-31.8 

-100.2 


UMAEAST 

40.7 

-18.1 

-88.1 

-216.8 

LOOP4 

39.9 

-22.9 

-100.2 

-204.7 

NGG3310 

30.2 

-30.2 

-76.0 

-150.0 

MRK9 

35.1 

-33.4 

-70.4 


AG 

47.9 

-31.8 

-72.8 

-153.3 

SUBA 

50.2 

-12.7 

-79.7 

-259.8 

MBM23 

50.4 

-13.3 

-72.8 


091346A 

39.8 

-23.0 

-90.0 

-140.4 

MRK421 

35.1 

-44.6 

-120.3 


GDFS 

41.9 

-19.9 

-60.1 

-119.8 


tainty of the emissivities. We note that this is not a test 
of uniqueness, but of consistency. But additionally, if the 
two emissivities were the same, then the motivation for 
finding a precise separation of LVC and I VC would be 
moot. 

For many emissivity analyses (e.g., Planck Collabo¬ 
ration XXIV 2011), some of the spectral data end up 
masked - mostly due to the presence of unaccounted hy¬ 
drogen in the form of H 2 or sometimes H+. In such cases, 
it might be beneficial to iterate on the standard devia¬ 
tion spectrum to determine the ideal separation between 
components including only the relevant data, those in 
the retained (not-masked) region. This has not been im¬ 
plemented here. 

Another consideration for closely adjacent ranges is 
that some signal from the LVC gas can contaminate the 
I VC component and vice versa, because of the overlap 
of the extended wings of their respective line profiles. 
An alternative method that takes into account these “in¬ 
truding” extended wings implicitly is based on Gaussian 
decomposition of the individual spectra (Section 8). 

6 . MAPS OF THE LINE INTEGRAL 

It is informative to calculate the integral Wm of the 
HI emission spectrum over velocity (or over distinct ve¬ 
locity ranges - e.g., LVC, IVC, HVC), because Hhi is 
related to the column density Vhi- At high latitudes, 
where the brightness of the infrared dust emission - the 
“cirrus” - is low because of low dust column density, the 


high-latitude 21-cm emission is also faint, with the peak 
temperature of the spectral lines usually small compared 
to the likely spin temperature Tg, so that the emission 
is optically thin. In this limit, the column density Vhi 
derived directly from IFhi is that corresponding to an 
infinite spin temperature: 

iVHi(oo)/C = Jndv = Whi , (1) 

where the conversion factor C = 1.823 x 10^^ cm“^(K km 
s“^)“^ and the integral is over the velocity range specific 
to a given component. This can be carried out for each 
spectrum in the data cube, producing a column density 
map for each of the chosen components. For optically 
thin emission, this is of course accomplished equivalently 
by summing channel maps. 

We have made Wm maps for each GHIGLS field using 
the velocity cuts in Table 2. Representative maps of Wm 
are shown in Figure 5 for the three velocity components 
in NEP. 

6.1. Nm 

The quantity Wm is a direct observable. Gomputation 
of maps of column density Nm with allowance for the 
effects of optical depth is discussed in Appendix G. For 
the low column densities characteristic of most GHIGLS 
fields, the corrections are small and so the Whi maps 
scaled by C are very close to those shown as Nm 
in Planck Gollaboration XXIV (2011) for many of the 
GHIGLS fields, aside from cropping, gridding differences 
(Section 3.2), and the effect the latter has on the base¬ 
line removal. In the results below and on the GHIGLS 
archive Nm has been calculated assuming Tg = 80 K 
(Equation (C2)). On the archive we have made avail¬ 
able EITS files of Nm maps for each field. Each file is 
a cube with five planes (0 to 4) corresponding to the ve¬ 
locity components HVG, IVG, and LVC and in addition 
IVC+LVC and HVC-j-IVC-l-LVC, respectively. We note 
that all of these maps are in units 10^^ cm“^. 

Eor each component map of the GHIGLS fields. Ta¬ 
ble 3 gives the mean Vri, the lower and upper ends of 
the range of Nm as defined by the 0.1 % and 99.9 % per¬ 
centiles, and the uncertainty ctnhi- 

Our estimates of ctnhi in the first of two columns are 
derived from differences in repeat observations. Where 
only single observations are available, ctnhi is the sum¬ 
mation in quadrature of noise, baseline, stray radiation, 
and scaling uncertainties (Boothroyd et al. 2011). In 
Appendix D we present a complementary assessment of 
the uncertainties in Vri, evaluating a second estimate 
of ctnhi with a power spectrum analysis (Appendix D.2). 
These values, in the second of the two columns, are in 
good agreement with those from the first approach. Un¬ 
certainties that could arise from correction for opacity of 
the line are discussed in Appendix C. 

In the tabulated units of 10^^ cm“^ the mean col¬ 
umn density in LVC ranges from 3.2 in UMl to 62.3 
in HS0624; in IVC, from 0.8 in MRK290 to 14.4 in NEP; 
and in HVC, from undetectable to 0.02 in LOOP4, to 0.1 
in UM2M, and to 6.2 in MG. Summing these, the mean 
total column density Nm ranges from 5.6 in UM2M to 
67.4 in HS0624. 

Even for the low column densities typical of the 
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Table 3 

Characteristics of Nm Component Maps (in 10^^ cm“^) 


Name 

(A^hi) 

Nm 

LVC 
Range ^ 

(t-Nhi ^ 

{Nm) 

A^m 

IVC 

Range 

o-Nhi 

(Nm) 

A^m 

HVC 

Range 

(t-Nhi 

MC 

6.6 

3.8 

11.4 

0.2 

0.1 

7.1 

4.4 

13.2 

0.1 

0.1 

6.2 

0.0 

21.2 

0.4 

0.3 

BOOTES 

7.0 

5.2 

11.4 

0.2 

0.1 

3.9 

1.3 

12.1 

0.2 

0.1 

0.2 

0.0 

1.9 

0.1 

0.1 

Necklace 

4.5 

3.7 

5.6 

0.1 

0.1 

4.4 

2.6 

5.9 

0.2 

0.1 






0X3 

10.9 

7.4 

15.9 

0.4 

0.4 

1.6 

0.0 

4.4 

0.4 

0.5 

1.6 

0.0 

7.6 

0.4 

0.7 

N1 

6.3 

4.2 

14.3 

0.2 

0.2 

3.1 

1.5 

15.1 

0.2 

0.1 

3.2 

0.0 

14.2 

0.2 

0.2 

G86 

11.0 

4.5 

25.1 

0.1 

0.1 

8.0 

1.3 

37.3 

0.2 

0.1 

0.4 

0.0 

2.9 

0.1 

0.1 

ELS 

17.6 

11.9 

26.9 

0.4 

0.2 

4.2 

2.0 

9.2 

0.4 

0.3 

2.0 

0.0 

7.5 

0.4 

0.3 

MRK290 

10.0 

7.0 

18.8 

0.4 

0.3 

0.8 

0.0 

9.6 

0.2 

0.1 

4.1 

0.0 

14.2 

0.5 

0.3 

DRACO 

6.5 

4.6 

16.3 

0.1 

0.1 

11.9 

4.2 

35.4 

0.2 

0.2 

3.7 

0.2 

16.0 

0.4 

0.3 

GROTH 

6.5 

4.8 

13.4 

0.3 

0.2 

3.1 

0.8 

4.7 

0.2 

0.1 

0.6 

0.0 

3.3 

0.3 

0.2 

NEP 

25.9 

12.5 

56.6 

0.3 

0.3 

14.4 

2.7 

56.7 

0.2 

0.2 

1.4 

0.0 

9.3 

0.4 

0.3 

UMIN 

52.0 

33.4 

67.9 

0.5 

0.4 

6.1 

3.5 

13.0 

0.2 

0.2 

2.3 

0.2 

11.7 

0.4 

0.2 

POL 

56.5 

15.8 

97.3 

0.9 

0.5 

10.3 

4.1 

21.0 

0.5 

0.3 

0.4 

0.0 

5.0 

0.6 

0.3 

POLNOR 

35.1 

8.8 

84.3 

0.4 

0.4 

7.8 

1.8 

16.9 

0.4 

0.3 

0.2 

0.0 

8.0 

0.6 

0.2 

MRK205 

23.4 

9.3 

57.2 

0.4 

0.3 

5.0 

1.6 

13.4 

0.4 

0.2 

0.2 

0.0 

9.9 

0.6 

0.2 

DEN 

4.4 

2.0 

6.5 

0.4 

0.2 

6.9 

4.0 

10.0 

0.4 

0.3 

2.6 

0.0 

9.7 

0.3 

0.3 

SP 

6.2 

3.8 

12.6 

0.2 

0.2 

3.7 

1.1 

9.6 

0.2 

0.2 

1.9 

0.0 

9.3 

0.3 

0.2 

UMl 

3.2 

1.0 

6.1 

0.5 

0.4 

13.2 

3.8 

38.4 

0.4 

0.4 

1.8 

0.0 

12.9 

0.4 

0.3 

SPIDER 

19.9 

6.2 

69.9 

0.3 

0.3 

7.3 

2.1 

16.9 

0.3 

0.3 

0.2 

0.0 

4.2 

0.2 

0.2 

SPC 

28.0 

9.9 

83.7 

0.5 

0.3 

2.6 

0.2 

10.9 

0.2 

0.2 

0.8 

0.0 

25.9 

0.4 

0.2 

1H0717 

32.4 

21.0 

70.8 

0.5 

0.4 

3.5 

0.9 

18.6 

0.4 

0.2 

0.7 

0.0 

12.1 

0.5 

0.2 

UMA 

27.5 

10.2 

74.4 

0.4 

0.3 

9.4 

3.7 

24.2 

0.4 

0.2 

0.9 

0.0 

31.6 

0.6 

0.2 

HS0624 

62.3 

31.5 

109.0 

0.5 

0.8 

2.6 

0.7 

6.6 

0.2 

0.2 

2.5 

0.0 

9.3 

0.4 

0.3 

UM3 

24.0 

4.5 

59.7 

0.4 

0.3 

6.6 

1.7 

19.7 

0.3 

0.3 

0.4 

0.0 

6.0 

0.4 

0.2 

09A079 

4.1 

2.4 

7.8 

0.5 

0.2 

1.6 

0.6 

3.0 

0.3 

0.1 

0.3 

0.0 

2.0 

0.4 

0.2 

UM2M 

3.3 

1.4 

7.2 

0.4 

0.3 

2.2 

0.4 

6.2 

0.3 

0.2 

0.1 

0.0 

1.9 

0.3 

0.1 

MS0700 

49.8 

30.0 

76.5 

0.5 

0.4 

2.3 

0.7 

6.5 

0.4 

0.2 






UMAEAST 

31.4 

17.7 

55.3 

0.4 

0.3 

9.7 

2.7 

23.4 

0.4 

0.3 

2.9 

0.0 

31.2 

0.7 

0.4 

LOOP4 

34.5 

22.7 

47.7 

0.4 

0.4 

8.7 

2.0 

17.7 

0.4 

0.3 

0.02 

0.0 

3.5 

0.6 

0.2 

NGC3310 

5.8 

2.1 

16.8 

0.3 

0.2 

4.1 

1.5 

7.7 

0.3 

0.2 

0.5 

0.0 

2.4 

0.4 

0.2 

MRK9 

43.9 

25.1 

67.6 

0.4 

0.3 

2.5 

0.6 

9.0 

0.2 

0.1 






AG 

5.2 

2.6 

11.8 

0.3 

0.2 

9.5 

4.3 

16.7 

0.2 

0.1 

3.9 

0.0 

28.4 

0.2 

0.2 

SUBA 

13.9 

11.4 

17.2 

0.3 

0.3 

6.8 

3.9 

11.1 

0.4 

0.3 

1.2 

0.0 

8.9 

0.9 

0.6 

MBM23 

57.0 

31.1 

89.8 

0.5 

0.5 

8.2 

3.7 

18.6 

0.3 

0.2 






091346A 

12.9 

9.3 

17.7 

0.3 

0.3 

1.9 

0.7 

4.7 

0.4 

0.3 

0.2 

0.0 

1.5 

0.3 

0.1 

MRK421 

8.2 

4.7 

14.4 

0.4 

0.3 

6.7 

2.4 

23.0 

0.4 

0.2 






CDFS 

5.6 

4.0 

9.7 

0.3 

0.4 

0.8 

0.2 

2.8 

0.2 

0.2 

0.2 

0.0 

3.5 

0.3 

0.3 


^ The two columns give the lower and upper ends of the range of Nm as defined by the 0.1 % and 99.9 % percentiles of values in 
the map. Negative Yhi, which shows up at a level equivalent to the Nm noise in faint HVC fields, is excluded in this calculation. 
^ First column: The uncertainty is calculated from repeat observations where available (Table 1), or from estimates of 

noise in the line emission, baseline, stray radiation, and scaling uncertainties following Boothroyd et al. (2011). Second column: 
Complementary calculation from power spectrum analysis (Appendix D.2). 


GHIGLS fields the data are of high sensitivity, with noise 
much less than the range of A^hi within the component 
maps. Those with the highest range of A^hi are the most 
useful for correlating with dust maps, because the mor¬ 
phological match will be better defined statistically. 

To identify possible issues in data calibration and re¬ 
duction, GET ACS data in a few of the targeted GHIGLS 
fields are compared to data from a new generation of 
wide-area H I surveys, in the north EEHIS (Appendix E) 
and in the south GASS (Appendix E). As shown by 
Eoothroyd et al. (2011), GHIGLS data agree well with 
the LAE H I survey data, in scale to within a few per¬ 
cent. We find that the agreement with EEHIS data is 
equally good. A more limited comparison with GASS 
reveals a calibration difference of about 6%. We note 
that comparison of GET 21-cm values of A^hi with those 
derived from measurement of Lya absorption toward 
high-latitude quasars are in good agreement, in the ratio 
1.00 ± 0.11 (Wakker et al. 2011). 

7. ANGULAR POWER SPECTRUM OF Nm 

The two-dimensional angular power spectrum of a map 
(image of the sky) f{x,y) is the square of the modulus 


of the Eourier transform f{kx,ky) where k is the spatial 
frequency (wavenumber) in the Eourier plane (Miville- 
Deschenes et al. 2007): 

P{k:,,ky) = \f{k:,,ky)f . (2) 

We apodize the image /(x, y), from which the median is 
removed, with a cosine function along its boundary prior 
to the Eourier transform operation in order to reduce 
edge effects, which otherwise produce a centered cross 
in the P{kx^ky) image (Miville-Deschenes et al. 2002). 
In practice, apodizing over five pixels at each edge is 
normally sufficient. 

The collapsed one-dimensional power spectrum P{k) is 
the azimuthal average of P(kx^ky) on a series of annuli 
of constant k = {k^ + Uncertainties are assigned 

based on the standard deviation of the mean within each 
of these annuli. If the above-mentioned cross is problem¬ 
atic, then adopting the median rather than the average 
is an effective alternative; in fact to be conservative we 
always adopted the median. To mitigate against these ef¬ 
fects further, the data point at the lowest k was excluded 
in the analysis below. 

The dots in Eigure 6 illustrate the basic anatomy of 
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Figure 5. Wm maps for NEP (upper, middle, lower: LVC, IVC, 
HVC). Note that the range in the colorbar in the HVC map is a 
factor of ten smaller so that the noise level is apparent in the map. 


the power spectrum of our data. At small k (large spa¬ 
tial scales in the map), P{k) roughly follows a power law: 
Po (k/ko)^ ^ where Pq is the amplitude of the power spec¬ 
trum at some representative scale ko and 7 is the scaling 
exponent. The exponent is alternatively called the spec¬ 
tral index or the slope (in a log—log representation as in 
Figure 6 ). The size of the exponent and its variations 
from one type of map to another provide insight into the 
turbulent structure of the ISM (Hennebelle & Falgarone 



0.001 0.010 0.100 

k (arcmin') 


Figure 6. Median power spectra of component maps of Nm/C in 
the NEP (black: LVC; green: IVC; and blue: HVC). Overlaid are 
fits (solid lines) to the model in Equation (5). As shown only for 
the case of LVC, this model consists of a power law (dash dot line), 
modified by the CBT beam (resulting in dashed line), plus the 
scaled noise template (dotted line). Shaded regions define ranges 
of k excluded from the fit. 


2012 ). 

At large k the power law is modified by the effect of 
the point spread function (beam) of the telescope and 
the noise. For a symmetrical Gaussian beam (j) described 
in one spatial dimension by its FWHM we have 

(l){x) = exp and 4>{k) = exp > (3) 

where the dispersion a = FWHM/(2^2 In 2) is related 
to <j/e by the Fourier relation (Jk = l/(27rcr). Above 
k ^ 0.04, where spatial scales are comparable to the 
size of the GET beam, the intrinsic P{k) is reduced mul- 
tiplicatively by <j)‘^{k). This can be seen in the power 
spectra in Figure 6 . At even higher k, the noise domi¬ 
nates; the distinctive shape of the noise power spectrum 
is discussed in Appendix D. 2 . 

We have investigated the effect of the asymmetrical 
beam. For the range of spatial scales sampled before the 
signal disappears into the noise, the optimal beam size 
is 9.'24 when treated as a Gaussian. As an alternative, 
starting with an image of the effective beam calculated 
for our modified Bessel function gridded data cubes (Sec¬ 
tion 3.2), we computed from its power spectrum the az- 
imuthally averaged 4)‘^{k). From our model of the GET 
beam we estimate an additional uncertainty 

SP{k) = bP{k){l-^'^{k)), (4) 

where b is an adjustable fractional error that we set to 
0.07. This is added in quadrature to the uncertainty de¬ 
rived for each annulus of constant k. While it is only 
an approximation, 6P has the desired effect of assign¬ 
ing lower weight to data strongly affected by the beam 
and still fitting the noise adequately. The parameters de¬ 
rived from the fits below are not sensitive to the precise 
choice of b. We also mitigated against further uncertainty 
by excluding data above a value /Cmax; again the results 
below are not sensitive to the precise value so long as 
the noise can be adequately assessed, and we adopted 
0.12 arcrnin”^. 
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7.1. Power Spectrum Model 

The computed power spectrum is fit to data iu the re¬ 
stricted rauge of k with a parameterized model cousistiug 
of a power law, modified by the beam, plus uoise: 

Pmodelik) = ^^{k) Po {k/koV + pN{k) . (5) 

The origiu of the uoise model, N{k)^ scaled by fittiug fac¬ 
tor 77 , is discussed iu Appeudix D. 2 . The best-fit param¬ 
eters are fouud usiug the above uucertaiuties as weights 
iu the IDL routiue mpfit.pro (Markwardt 2009). The 
derived model is plotted iu Figure 6 . We also iuspected 
a plot of P{k)/Pyj^o^Q\{k) to verify that there was uo bias 
hiddeu by the logarithmic display. 

For the three compoueut maps of A^hi of the NEP field 
the model fit to the data (see Figure 6 ) yields expoueuts 
- 2.86 ± 0.04, -2.69 ± 0.04, aud -2.59 ± 0.07 for LVC, 
IVC, aud HVC, respectively. The uucertaiuties cited are 
the formal Icr errors from the fits aud do uot iuclude auy 
systematic uucertaiuties. From exteusive testiug of al¬ 
ter uative choices iu the fittiug aualysis, we estimate that 
systematic uucertaiuties of the expoueut are of order 0.1 
for the LVC map aud probably somewhat smaller for the 
IVC aud HVC maps. We also fouud that relative dif- 
fereuces betweeu the expoueuts for differeut compoueuts 
are robust agaiust the systematic effects. The IVC com¬ 
poueut has a margiually shallower spectrum compared 
to the LVC compoueut, as cau be seeu directly from the 
data iu Figure 6 . 

With respect to the summary by Heuuebelle & Fal- 
garoue ( 2012 ), the expoueut for NEP LVC is roughly 
cousisteut with the value of —2.75 for H I seeu iu ab- 
sorptiou (Deshpaude et al. 2000), but by coutrast much 
shallower thau the expoueut of —3.6 for HI iu emissiou 
fouud by Miville-Descheues et al. (2003a) iu au iuterme- 
diate latitude field iu Ursa Major. The comparisou of 
power spectra compoueut by compoueut amoug differ¬ 
eut GHIGLS fields is au iuterestiug topic takeu up iu a 
forthcomiug paper. 

The values from H I columu deusity cau be compared 
to those from dust emissiou maps, for example the value 
of —2.9 for dust calculated for four regious at 100 /im 
(Gautier et al. 1992). Exteudiug aualysis to the eutire 
high latitude sky Miville-Descheues et al. (2007) fouud a 
mediau —2.93 with cousiderable dispersiou (0.3) from re- 
giou to regiou aud a treud of flatteuiug for faiuter regious. 
A value of —2.7 ± 0.1 was obtaiued for dust mapped 
iu the Polaris flare with Herschel/SPIME 250 /im aud 
IRIS 100 /im (Miville-Descheues et al. 2010). The Planck 
857 GHz mask-differeuced power spectrum (Plauck Col- 
laboratiou XVI 2014), characteriziug Galactic dust iu the 
iutermediate latitude sky, has au expoueut reachiug — 2.6 
asymptotically at high i multipoles, iu the rauge relevaut 
to our aualysis. 

We uote that the expoueut for the Planck 353 GHz EE 
aud BE polarized dust power spectra iu the iutermediate 
latitude sky (Plauck Collaboratiou lut. XXX 2014) is 
quite similar, —2.42, eveu though these spectra depeud 
additioually ou the liukage of the statistical properties of 
deusity aud maguetic field geometry. 

A detailed discussiou of the relatiouship betweeu the 
power spectra of Vhi maps aud dust emissiou maps of 
the same GHIGLS fields will also be preseuted iu the 
forthcomiug paper. 


8 . GAUSSIAN DECOMPOSITION OF LINE PROFILES 

With a goal of uuderstaudiug the physical couditious 
aud velocity structure iu the diffuse regious where these 
spectra origiuate, we have aualysed the H I hue profiles 
by decompositiou iuto a series of Gaussiau fuuctious us¬ 
iug a method similar to that of Hand (2000), details of 
which cau be fouud iu auother forthcomiug paper (Bla- 
grave, K. et al. 2015, iu preparatiou). Each spectrum is 
fit withiu the uoise with multiple Gaussiau compoueuts, 
with Gaussiau parameters amplitude, ceutroid velocity, 
aud EWHM (or dispersiou a). Each spectrum is fit iudi- 
vidually iudepeudeut of its ueighbors, but further itera¬ 
tive modificatious could iuclude iuformatiou from ueigh- 
boriug solutious, resultiug iu a more spatially cohereut 
set of compoueuts. 

As described further below (see also Blagrave, K. et al. 
2015, iu preparatiou), usiug various simulatious based ou 
Saury et al. (2014) we have fouud that the decompositiou 
distiuguishes effectively differeut parcels of gas at differ¬ 
eut temperatures. Erom simulated IVC plus LVC cubes 
(as iu Sectiou 5), we have also examiued various sets of 
two-dimeusioual histograms of the Gaussiau parameters 
iu order to differeutiate IVC aud LVC gas. 

8 . 1 . Distributions of Recovered Gaussian Parameters 

We illustrate the results ou GHIGLS data usiug the 
NEP field for which the average uumber of Gaussiaus 
fit to a giveu spectrum is 4.8. Eigure 7 shows a pair 
of two-dimeusioual histograms of the parameters for 
all Gaussiau compoueuts iu the EWHM—ceutroid aud 
EWHM—amplitude plaues. The EWHM-ceutroid his¬ 
togram (left) typically reveals what cau be called “pil¬ 
lars,” vertically-aligued features with au approximately 
coustaut ceutroid over a rauge of EWHM values. A 
particularly distiuct example of such a pillar appears at 
T ^ 0 kms“^ iu the rauge of the NEP LVC gas. Iu other 
GHIGLS fields, pillars cau be seeu at other LVC aud IVC 
ceutroid velocities, with varyiug degrees of coutrast. 

There is a sharp decrease iu the uumber of hues as the 
recovered EWHM decreases below 2 kms“^ (e.g., Eig¬ 
ure 7, right), similar to what has beeu fouud for both 
the LAB (Hand & Kalberla 2007) aud GASS (Kalberla 
& Hand 2015) surveys: H I emissiou hues uarrower 
thau EWHM 2 kms“^ are rarely fouud. A EWHM of 
2 kms“^ could arise as the thermal hue width of 100 K 
gas or might reflect the coutributiou of turbuleuce if the 
gas is eveu colder. The few compoueuts with EWHM 
< 1 kms“^ are uuphysical, the result of the automated 
fittiug routiue attemptiug to improve the model by fit¬ 
tiug Gaussiaus to rare remaiuiug uoise spikes. 

There are a few other features iu three-dimeusioual 
parameter space that we have ideutified as uuphysical 
Gaussiau compoueuts that arise because of the limita- 
tious of automated uucoustraiued Gaussiau fittiug rou- 
tiues. Eor example, a uumber of features with very low 
Hhi improve the model by fittiug Gaussiaus to resid¬ 
ual baseliues. These cau be excluded iu Eigure 7 right; 
auy compoueut with Whi < 1.0 K kms“^ is below the 
detectiou limit as determiued by the error aualysis iu 
Boothroyd et al. (2011) aud uucertaiuties iu Table 3. 

At the other extreme, there are very broad compo¬ 
ueuts with low amplitude that show up cousisteutly with 
ceutroids betweeu LVC aud IVC iu the EWHM-ceutroid 
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Figure 7. Histograms with respect to parameters of Gaussian components fit to spectra in the NEP field. The number of components 
in each two-parameter bin, /, is shown on a logarithmic scale quantified by the colorbar. Left: FWHM—centroid histogram. Components 
of the LVC gas with centroids near 0 kms“^ have a large range of FWHM and the ridge in the histogram forms a distinctive vertical 
“pillar.” At large FWHM the ridge in the histogram “arches” toward lower centroid velocities because of components with unphysically 
large FWHM fit to low amplitude emission from both IVC and LVC gas (see text). Right: FWHM—amplitude histogram. Lines of 
constant Wm (K kms“^) are shown as solid and dotted diagonal lines. Components with FWHM < 1 kms“^ (dashed line) are flagged 
as unphysical; these typically have low Wm. Any component with Wm < 1.0 K kms“^ is below the noise limit. 


histograms (Figure 7 left), creating “arches” between the 
pillars. Similarly, arches of extremely broad components 
also appear in the FWHM—centroid histogram with cen¬ 
troids between LVC and HVC. We have demonstrated 
that these too are unphysical, as follows. As in Sec¬ 
tion 5, we used the hydro dynamical simulations of Saury 
et al. (2014) to create an IVC plus LVC cube. Although 
this cube contains no intrinsic very broad components, 
we found on Gaussian decomposition (Blagrave, K. et 
al. 2015, in preparation) that such very broad compo¬ 
nents nevertheless do arise, again appearing at interme¬ 
diate velocities between the LVC and IVC pillars in the 
FWHM—centroid histogram. In the simulation and ac¬ 
tual data these are artifacts of the Gaussian fitting rou¬ 
tine, a result of fitting a single Gaussian to a merger of 
two (or more) weaker, narrower Gaussian components. 

8.2. Maps of Nui for Components with Different Line 

Widths 

The aforementioned LVG, IVG, and HVG velocity- 
selected components would be closely related to the 
Gaussian components forming vertical pillars within 
a range of centroid velocity in FWHM—centroid his¬ 
tograms. For example we would select Gaussian compo¬ 
nents as LVG if their centroids fall in the velocity range 
for LVG from Table 2, which for NEP is —20.5 < v < 
+47.9 kms“^. 

These histograms suggest a possible new complemen¬ 
tary direction for the subdivision of a cube into compo¬ 
nents on the basis of distinctive FWHM. As mentioned 
above narrow and broad components are commonly asso¬ 
ciated with two phases of the diffuse neutral atomic ISM: 
the GNM and the WNM, respectively. Analysis of his¬ 
tograms from all of the GHIGLS Gaussian components 
suggests that the natural division FWHM^ between nar¬ 
row and broad Gaussian components for the LVG emis¬ 
sion occurs at about 7.1 kms“^. With this division we 
can make Vhi maps corresponding to GNM and WNM 
by summing up the contributions from the appropriate 


components. To allow for uncertainties introduced by 
the Gaussian decomposition model and particularly by 
a sharp divide at FWHM^^ = 7.1 kms“^, we have em¬ 
ployed a Monte Garlo approach to produce and average 
~ 100 versions of the map using log FWHM^^ drawn from 
a normal distribution with a Icr dispersion of 0.05 dex 
about the indicated mean. 

The GNM-WNM separation using the Gaussian com¬ 
ponent approach can be tested using simulations in which 
the true results are known a priori because the gas tem¬ 
perature is known. For example, we have studied (Bla¬ 
grave, K. et al. 2015, in preparation) a simulation from 
Saury et al. (2014) for which ( = 0.2, Vs = 12.5 kms“^, 
and no = 1-0 cm“^, resulting in /cnm ^0.29. As we 
did above for the GHIGLS observations, in that work 
we selected a division in FWHM based on the distri¬ 
bution of Gaussian components for the simulation, in 
that case finding FWHM^^ = 4.5 kms“^. This results 
iR /cnm ^ 0.32, consistent with the actual /cnm for the 
simulation. This test also showed reasonable agreement 
between narrow and broad-component column density 
maps created from the Gaussian components and corre¬ 
sponding maps created from gas in known temperature 
ranges. 

The GNM Vhi map thus made from the selection of 
Gaussian components for NEP LVG is shown in Eigure 8, 
upper, as Nm/C. The GNM map has a low column den¬ 
sity compared to the total LVG but also highlights the 
filamentary structure running roughly diagonally across 
the field. The fraction of GNM in the LVG of NEP by 
column density is /cnm ^ 0.08, reflecting the dearth of 
narrow Gaussian components identified along many lines 
of sight in the NEP. We note that the GNM fraction by 
mass could be different, depending on the relative dis¬ 
tances of the GNM and WNM gas. This value of /cnm 
can be compared to those obtained directly from absorp¬ 
tion line spectra. Heiles & Troland (2003) found a global 
ratio of GNM to total Vhi of 0.39 for the Arecibo sky. 
This is an upper limit for the mass fraction of GNM be- 
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Figure 8. Column density maps for NEP for gas in the LVC 
range in units Nm/C from Equation (C2) after Gaussian decom¬ 
position of the line profiles. Upper: Column density from narrow 
components (FWHM < 7.1 kms“^) and assuming Tg = 80 K. 
Note that the range in the colorbar is a factor of five smaller than 
for the map below. This map can be interpreted as from the CNM. 
Filamentary structure is seen with a roughly diagonal orientation 
as in individual channel maps (cf. Figure 3, upper right). Lower: 
Column density map obtained by subtracting the above CNM map 
from the map of the total LVC component in Figure 5, upper. This 
is interpreted as being from the WNM and in this field is very sim¬ 
ilar to the total. 

cause of the systematic difference in distance between the 
CNM and WNM. Likewise, Dickey et al. (2009) found a 
fraction /cnm ^ 0.15 — 0.20 for the outer Galaxy. The 
result for CNM in the LVC range in NEP is therefore 
lower than these global values. 

We note that the observed low /cnm in NEP can be 
described well with appropriate tuning of the parameters 
in the Saury et al. (2014) simulations. Eor example, some 
simulations with ( = 0.3, Vg ^ 15 kms“^, and no = 
1.0 cm“^ result in /cnm =0.11. 

Also shown in Eigure 8, lower, is the WNM Vhi map. 
Because of the complications of the unphysical arches at 
high EWHM, this has been produced simply by subtract¬ 
ing the CNM map from that of the total LVC emission 
(Eigure 5, upper). The WNM map is quite similar to the 
total LVC map, given the low /cnm- 

8.3. Power Spectra of Maps of Vri for Components 
with Different Line Widths 



Figure 9. Power spectrum of Nm/C of the NEP as in Figure 6, 
but for LVC emission separated into CNM (red) and WNM (pur¬ 
ple) components (Figure 8). Overlaid as gray circles is the power 
spectrum for the total LVC Nm/C from Figure 6. The Gaus¬ 
sian decomposition introduces an additional uncharacterized noise 
component in the Vhi maps near the pixel scale, reflected at high 
spatial frequencies in the power spectrum above the modelled beam 
turnover (compare the purple and gray points); see text. 

The power spectra of the Vri maps of the LVC CNM 
and WNM components for NEP are shown in Eigure 9. 
The WNM component accounts for most of the Vri of 
the LVC emission and so its power spectrum is similar 
to that of the total LVC emission in Eigure 6. The CNM 
spectrum has less power overall and a shallower depen¬ 
dence on k. 

Note that the shape of these power spectra at high k 
(Eigure 9) is markedly different from that of the inte¬ 
grated Wri maps (Eigure 6). Despite attempts that can 
be made to keep solutions smoothly varying from pixel 
to pixel, the Gaussian fitting introduces an additional 
uncharacterized noise contribution near the pixel scale, 
which is reflected in the component Vri maps and prop¬ 
agates to a larger noise in the power spectrum at high k 
above the modelled beam turnover of the signal. In the 
model fitting we allow for this by the simple device of in¬ 
creasing the beam uncertainty parameter in Equation 4 
to b = 0.3, which in effect gives lower weight to the many 
data values at the highest k where noise now dominates. 

The model fits to the LVC WNM and CNM compo¬ 
nent maps (see Eigure 8) yield exponents of —2.7 ±0.1 
and —1.9 ±0.1, respectively. The exponent for the WNM 
component map is close to the exponent for the total LVC 
emission, —2.86 ± 0.04 (Section 7.1), within the system¬ 
atic errors, as expected because this component contains 
most of the mass. However, the power spectrum for the 
CNM component map is clearly shallower than for the 
WNM, a differential result that is robust against sys¬ 
tematic effects of how the model is fit. The significantly 
flatter spectrum found for the CNM map quantifies what 
can be foreseen readily in the maps themselves: the CNM 
map has more high-contrast small angular scale filamen¬ 
tary structure than the WNM map. 

Such a shallow dependence in a narrow-line compo¬ 
nent is not unexpected. The spectrum of the density on 
two-dimensional slices through the simulated cube used 
for the test in Section 8.2 has an exponent of —1.3, a 
much shallower dependence than from the Kolmogorov 
exponent of —8/3. Thus the spectrum of the three- 
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dimensional density would be —2.3 and furthermore this 
would be the exponent of the column density map of the 
cube (Miville-Deschenes et al. 2003b). 

Saury et al. (2014) also expect from their simulations 
of thermally bistable gas that the CNM power spectrum 
will be shallower than for the warm gas. Indeed, power 
spectra of A^hi maps produced from these simulations 
by selecting on the temperature of the gas (Blagrave, K. 
et al. 2015, in preparation) have much shallower power 
spectra for the CNM —1.9) as compared to the WNM 
(- -2.9). 

In addition to small scale structure corresponding to 
enhanced concentration of H I in the CNM, there are also 
possible effects from cold neutral gas becoming molecular 
and thus leaving structure in the remaining H I, albeit 
probably with lower contrast. The molecular transition 
is not modeled in these particular simulations. 

The CNM angular power spectra described here for 
both the NEP LVC range and the simulated data are 
significantly shallower than anything noted previously for 
H I emission, H I absorption, CO line emission, or dust 
emission (see Figure 10 in Hennebelle & Falgarone 2012). 

9. RELATIONSHIP OF H I CNM STRUCTURE TO THE 
ORIENTATION OF THE MACNETIC FIELD 

In Section 4 we commented on the striking filamen¬ 
tary structures in the LVC channel maps of NEP that 
cross the region roughly diagonally (lower right to up¬ 
per left). See the channel map in Figure 3. These fea¬ 
tures are fairly narrow in line width (present over only 
a few adjacent channel maps, as is evident in the movie 
of the cube) and so imprint on the Vhi map made from 
the narrow Gaussian components (Figure 8, upper). As 
is illustrated in Figure 10, this filamentary H I struc¬ 
ture in the condensed CNM in NEP is aligned roughly 
parallel to the direction of the Galactic magnetic field 
projected on the plane of the sky as inferred from the 
Planck 353 GHz thermal dust polarization map (Planck 
Gollaboration Int. XIX 2015). 

This relative orientation is in accord with the recent 
finding that the magnetic field tends to be oriented par¬ 
allel to the elongation of filamentary dust structures in 
both the high latitude sky (Planck Gollaboration Int. 
XXXII 2014) and in nearby Gould Belt molecular clouds 
(Planck Gollaboration Int. XXXV 2015). Such a sys¬ 
tematic tendency in relative orientation is important for 
understanding the observed power in B-mode relative 
to E-mode dust polarization (Planck Gollaboration Int. 
XXXVIII 2015). 

It is also interesting to note that the Gygnus spur end of 
Loop III seen in synchrotron emission also crosses NEP 
on the same diagonal and that the polarization of the 
synchrotron emission indicates that the magnetic field is 
parallel to the loop (Planck Gollaboration XXV 2015). 
Thus in NEP the predominant orientation of the field 
revealed by thermal dust and synchrotron polarization is 
similar. 

The observed Planck polarization arises from dust not 
just in the GNM gas, but in the WNM gas within which 

Another way to generate such a shallow power-law depen¬ 
dence is with cold and supersonic gas, as seen in simulations of 
isothermal high Mach number flows (e.g., Kim & Ryu 2005); how¬ 
ever, these conditions seem less relevant to the general interstellar 
medium than those in the simulations of Saury et al. (2014). 
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Figure 10. View of the magnetic held measured by Planck over¬ 
laid on a map of Nm from selected LVC gas in NEP. The “drapery” 
pattern, produced using the line integral convolution (LIC, Cabral 
& Leedom 1993), indicates the orientation of magnetic held pro¬ 
jected on the plane of the sky, orthogonal to the observed dust 
polarization (Planck Collaboration Int. XIX 2015). The colors 
represent Nm integrated over seven channels at the peak of the 
LVC emission (—6.4 to —1.6 kms“^), highlighting the structure 
of CNM gas. This illustrates how the orientation of the CNM gas 
filaments tends to follow the magnetic held across the upper part 
of the image, among several interesting relationships between the 
gas and the held (see text). 

the GNM is embedded. Thus the field orientation high¬ 
lighted in the GNM structure is more pervasive. (In ad¬ 
dition, dust in the IVG contributes to the total emis¬ 
sion, with about half the emissivity and perhaps some 
polarization.) The relationship between gas and mag¬ 
netic field is quite interesting in its complexity. For ex¬ 
ample, the field appears to wrap around a depression in 
the gas emission at (/,6) = (93°, 28°) (see also the total 
LVG in Figure 5, upper), suggestive of a bubble. Below 
that the field appears to wrap around the edge of the 
enhanced gas emission. 

10. CONCLUSIONS 

In this paper we present data from the GBT H I survey, 
GHIGLS. This deep/sensitive survey covers targeted re¬ 
gions of the intermediate Galactic latitude sky including 
several with distinctive IVG and/or HVG features. 

The data have been calibrated and corrected for stray 
radiation following Boothroyd et al. (2011). Data for 
most regions has an rms noise < 100 mK in a 0.8 kms“^ 
channel. The effect of 2I-cm line opacity on the calcu¬ 
lated column density Vhi and various other small sources 
of uncertainty in the GHIGLS measurements of Vhi are 
discussed. 

GHIGLS data agree well with the LAB H I survey data, 
in scale to within a few percent. We find that the agree¬ 
ment with EBHIS data is equally good. A more limited 
comparison with GASS reveals a calibration difference of 
about 6 %. 

We divide the H I emission into components with dif¬ 
ferent velocities and produce Vhi maps of LVG, IVG, and 
HVG components. Structure in these maps is quantified 
by the angular power spectrum modelled with a power 
law. For the NEP field, the power-law exponent for LVG 
(—2.86 ± 0.04) is close to that found for dust maps in 
the intermediate Galactic latitude sky; the exponents for 
IVG and HVG reveal a marginally flatter power-law be- 
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havior (exponents —2.69 ±0.04 and —2.59 ±0.07, respec¬ 
tively) . 

We fit the spectral line profiles with multiple Gaussian 
components to differentiate gas characterized by different 
line widths, enabling a subdivision into emission by two 
ISM phases, the WNM and CNM. The CNM map of 
A^hi is unique in its angular power spectrum, having a 
power-law exponent of —1.9 ± 0.1 for LVC in the NEP 
field, shallower than seen in any component map but 
consistent with the power-law exponent of CNM seen in 
simulations. The flatter power-law behavior reflects more 
small scale structure associated with this phase. 

There is evidence that filamentary structure in the H I 
CNM gas is aligned with the Galactic magnetic field re¬ 
vealed by Planck polarization. 

GHIGLS data have already been used productively in 
a number of Galactic and extragalactic applications and 
should be interesting for many more. On the GHIGLS 
archive (www. cita.utoronto. ca/GHIGLS), fully-reduced 
data cubes, along with movies and Nm component maps, 
are available for inspection and downloading. 
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APPENDIX 

A. DATA OBTAINED USING THE GBT SPEGTRAL 
PROGESSOR 

As indicated by the notes in Table 1, 11 fields in¬ 
cluding the central portion of the BOOTES field were 
observed with the now-retired GBT Spectral Processor 
(GBT SP). As is the case for the GBT ACS data, these 
data were collected in in-band frequency-switched mode, 
except for the ELS field observed with out-of-band fre¬ 
quency switching. The ELS data were regridded from 
the observed 0.5 kms“^ spaced channels to the adopted 
ACS spacing of 0.807 kms“^, using a cubic spline inter¬ 
polation. Eor the other fields the channel spacing was 
maintained as observed at 1.03 kms“^ (for the entire 
BOOTES data, this was adopted even for the ACS por¬ 
tion). The number of repeats and the resulting emission- 
free channel noise, def, are found in Table 1. 

We processed the raw GBT SP data following the GBT 
ACS pipeline described in Section 3.1 and Boothroyd 
et al. (2011), including the benefit of the improved stray 
radiation correction over Lockman & Condon (2005). Eor 
consistency with the GBT ACS cubes, the spectra were 
gridded on a 3.'5 grid using the modified Bessel function. 

There is a slightly different calibration unique to the 
GBT SP because in that era there was a different noise 
diode in the GBT 21-cm receiver. Therefore, the calibra¬ 
tion discussed in Boothroyd et al. (2011), specifically the 
scaling of Ta by the factor 1.024 ± 0.009, is not relevant. 


The scaling of Ta appropriate to the GBT SP data was in¬ 
vestigated in two ways. The first approach used archival 
GBT SP data on the S6 and S8 calibration standards. 
Eor these, the XX and YY spectra are slightly inconsis¬ 
tent with each other - differing by 5 to 10 % - but their 
average spectrum, (XX ± TT)/2, is reproducible and 
indicates consistently that the GBT SP measure of the 
main beam temperature Tb is 1.05 to 1.08 larger than for 
the GBT ACS spectra of these calibrators. The second 
approach compared regions of overlap between spectral 
cubes made with GBT SP observations with cubes made 
with GBT ACS observations: the original SP central field 
in BOOTES with the ACS flanking fields, UM2M with 
NGC3310, UM2M with 09A079, UM3 with UMA (and 
LISZTA), and ELS with NEP (see Eigure 1). These com¬ 
parisons indicate overestimates of the same order (1.05 
to 1.08). Therefore, in the final pipeline processing the 
GBT SP Ta data were scaled using a factor 0.97 (± 0.01) 
instead of 1.024 to bring the GBT SP data onto the same 
Tb scale as the GBT ACS data. 

B. ARGHIVAL DATA OBTAINED FOR 09A079 USING THE 
GBT AUTO-GORRELATION SPEGTROMETER 

The low column density region within the Lockman 
Hole is covered in part by the H I survey UM2M using 
the SP (Appendix A) and an overlapping field 09A079 
observed by Grossan et al. (2012) using the ACS (pro¬ 
posal GBT/09A-079). The latter data were obtained to 
estimate the Galactic dust foreground for studies of the 
CIBA using a 160 jam map from Spitzer and to probe 
the relation between dust and gas in this very low dust 
regime. Because of the low column density it is essen¬ 
tial to correct the H I spectra for stray radiation be¬ 
fore embarking on the science analysis. Therefore, we 
have reprocessed these ACS data using the pipeline that 
we developed (Section 3). This cube is available on the 
GHIGLS archive and we have used it in the LHM mosaic 
described in Section 3.3. These observations were carried 
out by mapping several overlapping subfields, some of 
which involved scans along constant Galactic longitude 
rather than the more usual constant Galactic latitude. 
The subfield structure and irregular coverage is not de¬ 
tailed in Table 1 , but is imprinted in the noise and weight 
maps. 

C. OPTIGAL DEPTH EFFEGTS 

A diagnostic warning that an observed H I spectrum 
might be affected by opacity is if Tb approaches the (plau¬ 
sible) spin temperature Tg. In the extreme the line pro¬ 
files might become flat topped, but that is rarely encoun¬ 
tered here. We wish to explore more subtle effects. 

In the absence of absorption-line measurements to 
combine with the emission measurements to distinguish 
the effects of optical depth r (Strasser & Taylor 2004; 
Dickey et al. 2009), as is the case for GHIGLS data, 
it is nevertheless possible to make an estimate of the 
potential optical depth effects for different assumed 
Tg. As a reasonable spin temperature we adopt Tg = 
80 K which is the collisional temperature found from 
intermediate-latitude H 2 observations for column densi¬ 
ties near 10^^ cm“^ (Gillmon et al. 2006; Wakker 2006). 
This would be appropriate for CNM. The CNM equilib¬ 
rium temperature is not constant but depends inversely 
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on the local density. Estimates of the CNM tempera¬ 
ture in the Galactic plane do show lower values (Dickey 
et al. 2003); however, Heiles & Troland (2003) find that 
for \b\ > 10° and CNM column densities near 10^^cm“^, 
the mass weighted Tg is 70 K (median) and 108 K (aver- 
age). 

C.l. Complications in a Multi-phase ISM 

Differences among and within individual observed H I 
line profiles for different lines of sight and their overall 
complexity are reminders that a single Tg is naive. The 
emission is summed over gas in different stable phases 
of atomic gas with different temperatures, the CNM and 
WNM, and there could be substantial amounts of warm 
thermally unstable gas as well (Heiles & Troland 2003; 
Hand & Kalberla 2007; Hennebelle & Falgarone 2012; 
Saury et al. 2014). While Tg for the CNM might be near 
80 K, Tg for the WNM is much higher (> 300 K; Dickey 
et al. 2009). One could carry out component (phase) 
separation by line profile decomposition (see Section 8), 
to segment CNM and WNM and then assign different Tg 
in correcting these classes of profile. But even this would 
not capture the complexities arising from the relative ge¬ 
ometries of the gas in different phases. 

However, in GHIGLS the highest values of Tb in the 
observed profiles are caused by the CNM. We have shown 
this by making a cube of the CNM gas emission using the 
Gaussian decomposition method discussed in Section 8. 
For each voxel we then formed the ratio Tb(CNM)/Tb 
and found that the ratio approaches unity for higher Tb. 
Because optical depth corrections become important for 
Tb comparable to Tg, which is accentuated in CNM gas 
where Tg is low, such high-Tb parts of the profile are most 
in need of correction, using a Tg appropriate to CNM. 

On the other hand at intermediate latitudes the WNM 
produces broad profiles that have low Tb and so the mag¬ 
nitude of a correction even using a Tg that is inappropri¬ 
ately low will be inconsequential. Therefore, a simple 
correction using a single Tg appropriate to CNM might 
not be unreasonable. That is the approach that is eval¬ 
uated quantitatively here. 

C.2. Impact on Nm 

The actual column density is always larger than the 
direct line integral A^hi(oo) in Equation (1) and within 
our assumptions is given by 

Nn{T,)/C = jT,Tdv (Cl) 

= T,J-log{l-n/T,)dv (C2) 

= W^Hi + 1/(2T,) jT^dv + ... (C3) 

The first-order correction, l/(2Tg) f d^ from Equa¬ 
tion (C3), is quadratic in Tb. Therefore, the cube, 

and maps of integrals made from it, reveal where there is 
sensitivity to optical depth. The form of the first-order 
correction has the great utility of showing the explicit 
scaling with l/(2Tg), quantifying the importance of an 
appropriate choice of Tg. 

In Figure 11 we show the first-order correction from 
Equation (C3) as a function of the optically thin solu- 



Figure 11. Contours of histogram of the first-order correction, 
l/(2Ts) f di;, as a function of the optically thin solution VEhi 
for NEP, here using Tg = 80 K. Values for the three components 
are overplotted in different colors: LVC (black), I VC (red), and 
HVC (blue). A line segment of slope two reveals the approximate 
quadratic dependence. The unit-slope dotted lines indicate differ¬ 
ent percentage corrections VEhi- Horizontal colored dashed lines 
indicate the uncertainties ctwhi ~ ‘^Nhi/^ three compo¬ 

nents in the NEP field, from Table 3. 

tion Hhi pixel by pixel in NEP, distinguishing each of 
the three velocity range components. This log—log plot 
reveals an approximate quadratic dependence of the cor¬ 
rection in Wm (such a line is plotted), so that the frac¬ 
tional corrections scales as Whi, potentially rising above 
the uncertainty The spread in this figure is con¬ 

sistent with expectations from simple model spectra us¬ 
ing Gaussian components. The vertical height at a given 
Hhi is inversely related to the line width and the vertical 
spread is also influenced by the amount of overlap of the 
components in velocity. 

This correction can be compared to the typical uncer¬ 
tainty in A^hi discussed in Appendix D. For NEP a value 
cfnui of 0.2 X 10^^ cm“^ from Table 3 corresponds to Hhi 
^ 1 K kms“^. Any opacity correction below this limit 
can be ignored. Above this limit, LVC is affected at 
the 2—8 % level; the uncertainty in the correction will of 
course be less than this. Some I VC is affected in this 
same range but other I VC much less. The HVC com¬ 
ponent is least affected by any opacity (at the lower left 
the upturn relative to quadratic relates to the positive 
definite nature of the correction even for low Tb). 

In Figure 12 we present maps of this first-order cor¬ 
rection for the LVC and IVC components in NEP. The 
HVC correction, not shown, is minimal and below the 
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Figure 12. Map of the first order correction, l/(2Ts) f T^dv, for 
NEP (upper: LVC, lower: IVC). Here we set Tg = 80 K. Note that 
the range in the colorbar is smaller by a factor of ten compared 
to that for the corresponding Wui maps in Figure 5. The two 
contours show the (full) correction as a percentage of Wm (2.5 % 
and 5%). The HVC correction is minimal and so is not shown. 

noise (Figure 11). These appear like non-linear represen¬ 
tations of IThi in Figure 5. Also note the relative scales. 
A few contours show the (full) correction as a fraction 
of Wht Thus from Figures 11 and 12 both LVC and 
IVC components contain regions where opacity in the 
line is non-negligible and thus an uncertainty arises in 
the derivation of Vhi , dependent on the choice of Tg and 
more profoundly on the assumption of a constant Tg. 

C.3. Assessment using all GHIGLS Spectra 

Given the approximate quadratic dependence seen in 
Figure 11, it is advantageous to divide out one power of 
the dependence of the opacity correction on IFhi , so that 
the fractional correction (A^Hi(Fg)/C—1 Fhi)/1Fhi is seen 
to be approximately linear in Whi- This is plotted in 
Figure 13 for a single Tg = 80 K as a two-dimensional 
histogram using the values for every pixel in our com¬ 
ponent IFhi maps of all fields. The various islands, or 
ridges, within the two-dimensional histogram arise be¬ 
cause the line profiles leading to a given IFhi can be 
different: broad with low Tb and so a smaller correction, 
or narrow and peaked and so more affected. 

Also shown in the solid contours is the first-order cor¬ 
rection for the same Tg. The agreement is excellent at 
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Figure 13. Histogram of full fractional correction 

{Nm{Ts)/C-Wm)/Wm calculated for Tg = 80 K for all 

pixels in our component Hhi maps of all GHIGLS fields. Upper: 
LVG; middle: IVG; lower: HVG. A common line segment of 
slope one, transferred from Figure 11, reveals the approximate 
linear dependence. The solid contours, representing the first-order 
correction for the same Tg, show that this approximation remains 
quite good even at the highest values of Wm encountered. The 
dashed lines, also from Figure 11, are for typical values of the 
uncertainties in Nm/C for those components (from the 

NEP values in Table 3). The correction in the HVG component 
is rarely above this uncertainty. IVG and LVG are progressively 
more affected, having typically higher VEhu and so making an 
appropriate and accurate correction is more of a concern. 
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low Whi as expected and remains quite good even at 
the highest values of Wui in these intermediate latitude 
fields. This difference is much less than the uncertainty 
that arises from which value of Tg to adopt. 

We have repeated this, but for Tg = 200 K. To verify 
that the scaling is approximately inversely proportional 
to Tg for the full correction, we have multiplied the latter 
correction by 200/80 and overlaid contours on a figure 
like Figure 13. Again the agreement is excellent at low 
Whi as expected and it remains quite good even at the 
highest values of Whi encountered. 


C.4. Impact on Applications of Nm Maps 

As mentioned above, maps of A^hi have a variety of ap¬ 
plications. For example, in Planck Collaboration XXIV 
(2011) the GHIGLS Vhi maps were used in conjunction 
with maps of thermal dust emission to examine regional 
differences in the dust emissivity, the amount of emission 
per Vht To identify lines of sight along which there is 
a significant amount of hydrogen in molecular form, for 
which Vhi would no longer be representative of the to¬ 
tal column density, a masking procedure was developed. 
This was based empirically on the presence of emission 
in excess of that expected from the correlation at lower 
Vhi in a particular field. In retrospect the masked re¬ 
gions include all lines of sight in which GO is detected 
using the Type-3 Planck GO map (Planck Gollabora- 
tion XIII 2014), but the masked region is larger than 
this, which is interpreted as the presence of molecular, 
but GO-dark, gas for which the H I diagnostic signal is 
therefore “missing.” Qualitatively, there might be some 
ambiguity because Vhi could be underestimated if not 
corrected for opacity and indeed the masked regions in¬ 
clude areas with bright Tb and significant opacity such 
as in Figure 12. However, the Vhi maps used were al¬ 
ready corrected for opacity using Tg = 80 K. In principle, 
one could account for the excess dust emission by low¬ 
ering Tg below 80 K on a pixel-by-pixel basis, produc¬ 
ing a map of Tg in the masked region. In this solution, 
many lines of sight would require Tg to conspire to be 
as low as the peak Tb in the spectrum. Assuming that 
this is not allowed and that Tg = 80 K is a reasonable 
value for the entire field as justified earlier in this sec¬ 
tion, the excess dust emission cannot be accounted for 
by an insufficient opacity correction. Furthermore, the 
purpose of the masking was to exclude regions in which 
Vhi might not be a good tracer of total column density, 
and Vhi in the retained regions is certainly not sensitive 
to the choice of Tg. Thus by excluding regions with the 
masking procedure, the analysis of dust emissivity in the 
atomic gas is not significantly impacted by uncertainty 
in the opacity correction. 

However, often the full Vhi map is desirable, as in the 
power spectrum analysis. The appropriateness and im¬ 
pact of using the corrected Vhi in such a situation can 
be guided by a sensitivity analysis to gauge the robust¬ 
ness of the derived results when different values of Tg 
are adopted. For example, in the power spectra of Vhi 
corresponding to the analysis of Vhi component maps 
in NEP in Section 7.1 and Figure 6 with Tg = 80 K, if 
no correction is applied then the power spectra are only 
minimally steeper (< 0.02). 


D. ASSESSMENT OF UNCERTAINTY IN A^hi 

As discussed by Boothroyd et al. (2011), there are a 
number of distinct contributions that affect the accuracy 
of the spectra: noise, baseline fitting errors, and errors in 
the subtraction of the stray radiation spectrum. Because 
the column density is such an important data product for 
subsequent analyses, we focus here on the related ques¬ 
tion of the uncertainty ctnhi ^hi (for this discussion, 
Vhi is used interchangeably with IFhi, other than the 
scaling parameter C as in Equation (I)). The results 
have been given in Table 3. 

The following subsections summarize our understand¬ 
ing of the behavior of the various contributions. This 
exercise also involves various checks on the quality of the 
GHIGLS data, including a critical look for systematic 
errors. 

D.I. Noise in Emission-free Channels 

In emission-free channels, the noise def can be mea¬ 
sured as the rms fluctuations about zero in a spectrum 
or, alternatively, in a single channel map in the cube. 
The same noise can be assessed from the width of the 
peak near zero in a histogram of T\y{x,y,v). Eor our 
basic observation and processing - a single mapping of 
a field with 4 s integration per spectrum and modified 
Bessel interpolation into a cube with 0.8 kms“^ chan¬ 
nels - the measured noise, do, is typically 110 mK. With 
a total of Um repeats of the mapping, this is found to go 
down as cro/y/n^ as expected. Measured values of the 
typical (Jef for various GHIGLS fields are summarized in 
Table I. 

Because we are working with column density Vhi 
here, we are interested in how the noise increases when 
Uc channel maps are summed. Eor Uc independent 
emission-free channels, for a single mapping with do ^ 
no mK the expected rms error in Vhi would be 0.14 x 
IO^^^nc/75 cm“^ and for the GHIGLS observations this 
can be scaled by cref/cro(< I). This error is consistent 
with what is found directly from the rms of maps of the 
sum of emission-free channels (see also Section D.2.I). 
The choice of 75 channels for the illustrative normaliza¬ 
tion is a reasonable one, being close to the median num¬ 
ber of channels in the component intervals in Table 2; it 
represents a velocity interval of 60 kms“^. 

D.2. Power Spectrum of the Noise 

As is apparent from the dots in Eigure 6, the shape of 
the power spectrum of the noise is not flat (i.e., white) 
cA k > O.Iarcmin”^, but instead has a characteristic 
decrease at the largest k that arises from correlations in¬ 
duced by the modified Bessel function interpolation of 
each observed H I spectrum onto the grid (Section 3.2). 
We have verified the precise form of this decrease through 
simulations, starting with a white noise channel map that 
is subsequently spatially filtered (convolved) to corre¬ 
spond to the modified Bessel function interpolation used 
in the gridding of our data. 

When we sum a number of emission-free channels in 
our observed cube and compute its power spectrum, at 
high k the shape of this power spectrum is the same 
regardless of the number of channels. Eurthermore, this 
same shape is seen in the A^hi data itself (Eigure 6). 
Therefore the noise can be quantified simply by a “noise 
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Figure 14. Estimate of the noise in each channel of the NEP 
cube from the power-spectrum analysis, ^yrj = CTps/fJef (Ap¬ 
pendix D.2.1). At the ends of the spectrum ~ 1, as expected 
for primarily emission-free channels. However, where there is signal 
the estimated noise is somewhat larger. To model the contribution 
of noise arising from the 21-cm emission itself, we use the approxi¬ 
mation 1 + rmsTb(r’)/20 K, based on the rmsT^ measured in each 
channel map (Appendix D.3). The resulting smooth curve (right 
axis) captures some basic features of the increase in noise. 

template,” N{k), and a scale factor, r]. This is reflected 
in the full parameterized model of the one-dimensional 
power spectrum in Equation (5). We use the fitted level r] 
as a quantitative diagnostic of the noise (from all sources) 
in the map. We denote as dps the estimate of the rms 
noise in a map that is derived using this power-spectrum 
analysis. From this complementary approach estimates 
of ctathi for the component A^hi maps of the GHIGLS 
fields have been entered in Table 3. 

D.2.1. Noise Template 

In this paper the noise template used is that derived 
from emission-free end channels. This can be quantified 
by cTef, the rms noise of the inverse Fourier transform of 
this noise template. Because the same template, simply 
scaled by 77, is used in the power-spectrum analysis we 
have cTps = y/r^cTef. 

As an example, power-spectra of the NEP cube can 
be modelled on a channel by channel basis using Equa¬ 
tion (5) to find 77. The resulting spectrum = cr^s/o'ei 
is shown in Figure 14. This ^yrj spectrum approaches 
unity in the velocity ranges —350 < v < —250 kms“^ 
and 175 < v < 275 kms“^ because these are the ranges 
of the typical end channels used to remove baselines. The 
end channels used vary from spectrum to spectrum across 
the field. In these velocity ranges over 99 % of the pixels 
in each channel have been used to determine a baseline. 
This also demonstrates that the noise template is defined 
consistently. 

As another application, we made A^hi maps from an 
increasing number of channels of a simulated noise cube 
and found the expected behavior 77 oc Tic, or cTps oc y/n^. 
Fitting baselines has a notable effect on the accumulation 
of noise. To demonstrate this we fit baselines in this noise 
cube, using data where the emission-free end channels 
would normally be, and subtracted them. 

For A^hi maps made from channels where baselines 
were fit the noise 77 first increases with ric but eventuaiiy 
reaches a maximum and turns over when ric is compara- 


bie to the range used to fit the baselines, the exact point 
depending on the degree of the polynomial baseline being 
fit. This effect is a result of the correlations introduced 
by the baseline fitting. We examined the NEP cube in 
the same way and found a similar behavior with increas¬ 
ing Tic. 

For the channels where there is signal, away from where 
the baselines are defined, the noise is somewhat larger 
than simply def (Figure 14). This will affect the noise 
in A^hi maps too; the noise estimated from the sum of 
emission-free channels will always be an underestimate. 

D.3. Noise in the Line Emission 

The 21-cm line emission itself increases the noise by 
a factor approximately 1 + T\y{v)/T^y^ (Boothroyd et al. 
2011) and because the GBT 21-cm receiver noise is so 
low (Tgys ^ 20 K) the increase can be significant even in 
the GHIGLS fieids. 

This additional noise will not be uniform across the 
channel map, but the size of the effect can be demon¬ 
strated using a representative T\^{v). For each channel 
map we adopted the rms value, rmsTb(7;), noting that 
this spectrum is the sum in quadrature of the mean and 
standard deviation spectra seen in Figure 4. The ap¬ 
proximate increase in the noise spectrum predicted on 
this basis, 1 + rmsTb(7;)/20 K, is shown for NEP in Fig¬ 
ure 14. It captures basic features of the increase seen in 
the actual noise spectrum. 

D.4. Baseline Errors 

Fitting and subtracting a baseline introduces errors in 
A^hi through uncertainties in the fitted coefficients. Base¬ 
line errors add to the increased noise apparent in the 
shoulders adjacent to the emission-free channels in Fig¬ 
ure 14 and to the channels with more significant emission 
as well. 

This will propagate to uncertainty in A^hi maps too. 
Unlike the noise in the line emission, baseline errors 
are similar for adjacent channels in a spectrum and 
so accumulate as tt-c- Boothroyd et al. (2011) showed 
that these amount to an uncertainty in Nm typically 
0.3 X 10^^(71c/75) cm“^, comparable to the noise esti¬ 
mates above for tt-c ^ 75 but with a different Uc depen¬ 
dence. This was found by comparing independent XX 
and YY polarization cubes. These cubes have common 
stray radiation and so their difference (XX—YY) reflects 
the quality of baselines plus thermal noise (which can be 
seen increasing where there is a signal). 

The assumption made in both Planck Collaboration 
XXIV (2011) and Boothroyd et al. (2011) is that the 
noise in the difference and the noise in the sum are iden¬ 
tical. This is not the entire story, as we demonstrated by 
examination of the power spectra for a single observa¬ 
tion of a subfield of NEP. We produced maps of Wm 
from both XX-j-YY and XX—YY cubes. The power- 
spectrum analysis yielded slightly larger values of dps in 
the XX+YY case, by up to 25 %. 

D.5. Errors in the Stray Radiation Speetrum 

Boothroyd et al. (2011) showed that there is a sig¬ 
nificant error contribution to Vhi from the uncertainty 
in the predicted stray radiation spectrum that is sub¬ 
tracted. This is an error that can vary slowly over a field 







22 


Martin et al. 


in the course of mapping and like the baseline removal 
could have a non-linear effect on the scaling parameter 
of the noise power spectrum. 

An empirical approach to determining this uncertainty 
is to compare the reproducibility of A^hi derived from 
cubes made from different independent mappings of a 
field (up to three in our survey). This includes all sources 
of error but is dominated by the accumulated baseline 
errors and the imperfect stray radiation subtraction. In 
principle it would be possible to schedule observations to 
make repeated maps with virtually the same stray radi¬ 
ation. However, in practice because of our exploitation 
of flexible scheduling the mapping of (parts of) any field 
took place at different times of day, on different days, 
and at different times of year, and so generally the stray 
radiation and the estimated correction are different in 
different maps of the same field. Indeed, this was used 
by Boothroyd et al. (2011) to study the properties of the 
stray radiation and to calilDrate the amplitude of the side- 
lobes in the all-sky response of the GBT. A corollary is 
that because the stray radiation can now be predicted in 
advance, future observations could be scheduled for times 
at which the stray radiation is minimized. The data re¬ 
ported in this paper for the faint N1 field approach this 
ideal. 

For every GHIGLS field in addition to the observed 
emission-line cube we have a predicted stray radia¬ 
tion cube from which we can compute A^nistray Using 
the above difference approach, Boothroyd et al. (2011) 
showed that the error contribution to A^hi is of order 
0.07 A^HIstray 

Because the errors in two repeated observations i and 
j are independent, and ultimately the two are averaged 
together to form (A^h) = (^H,i + A^H,j)/2, the estimator 
of interest in assessing errors in A^hi is the dispersion 
(the standard deviation about the mean) of the map 
of AA^h = (^H,i — A^H,j)/2, or for the complementary 
approach being used here dps from the power-spectrum 
analysis of that difference map. 

Again using the power-spectrum analysis to find dps 
we are able to look at not only the noise in the repeat 
observation difference maps, but also the noise in their 
sum. The disagreement between the results is somewhat 
greater than in the analysis of XX and YY used to study 
baseline errors. The power-spectrum analysis yielded dps 
values for the sum roughly 50 % larger than for the dif¬ 
ference. 

E. COMPARISON WITH DATA FROM EBHIS 

A thorough comparison of our GHIGLS data with that 
from the LAB survey was performed in Boothroyd et al. 
(2011). Here we compare our data with that from the 
Effelsberg-Bonn H I Survey (EBHIS) (Elder et al. 2010; 
Winkel et al. 2010; Kerp et al. 2011; Winkel et al. 2012) 
with the goal of testing the quality of the GHIGLS data 
and the data reduction procedures, especially the cor¬ 
rection for stray radiation. The EBHIS collaboration 
has provided unpublished data cubes and stray radiation 
cubes for the NEP field that has been used for many il¬ 
lustrations in this paper, another large field centred on 
SPIDER but including portions of POL, SP, SPG, and 
UMA (these all relate to the NCP loop and so we use the 
designation NGPLEB here), and DRACO. 


E.l. EBHIS Data and GHIGLS Modifications 

The EBHIS data have an angular resolution (FWHM) 
^ 10.'5 (Winkel et al. 2010) on a Galactic SFL grid with 
3.'2564 pixels (we regridded to the GHIGLS 3.'5) and a 
channel spacing of Av = 1.288 kms“^. The noise maps 
for all three cubes are presented in the first column of 
Figure 15 where the EBHIS maps have been cropped to 
match the coverage of the equivalent GHIGLS field. The 
noise in each map pixel is set by the rms of emission- 
free channels in the corresponding spectrum. The mean 
noise is agf ^ 95 mK (NEP), ^ 86 mK (NGPLEB), and 
- 90 mK (DRACO). 

For a direct comparison with EBHIS, the GHIGLS 
data cubes first need to be convolved to the slightly lower 
resolution of the EBHIS survey. Following the convo¬ 
lution of each channel in the cubes, the GHIGLS data 
were regridded both spatially (cubic convolution interpo¬ 
lation using the IDL routine interpolate) and along the 
velocity axis (linear interpolation). These steps of con¬ 
volution and regridding reduce the original noise of the 
GHIGLS data for NEP, NGPLEB, and DRACO from 
68 mK, 75/105 mK,^^ and 61 mK (Table 1) by about 
a factor 2 to an average rms noise 34 mK, 40/55 mK, 
and 31 mK, respectively. The rms noise maps for the 
modified GHIGLS data are shown in the right column of 
Figure 15. 

The scanning and map-making strategies are reflected 
in the various geometric (grid) patterns (called GPs be¬ 
low) that appear in these noise maps. In addition, spec¬ 
tra with very few emission-free end channels due to the 
presence of galaxies result in noise peaks in these maps. 

E.2. Relative Galibration 

Prior to detailed comparisons between EBHIS and 
GHIGLS, we first determine if correction by a scale factor 
is necessary. This is done by comparing Tb along every 
common line of sight, fitting the scatter plots using the 
anticipated linear model: 

EBHIS = a X GHIGLS + b . (El) 

This is repeated for scatterplots of Wm as well. Both 
correlations are very good. 

The scale factor (slope) a should be close to unity, but 
reflects the different methods by which the spectra have 
been calibrated. In the ideal case, there would be no off¬ 
set b. However, inconsistencies remaining after the stray 
radiation and/or baseline corrections can introduce an 
offset and this can in turn influence the slope. These in¬ 
consistencies and noise will have the largest relative effect 
where Tb is small. To ensure that the derived scale factor 
is not influenced unduly we limit the data to Tb > Tnmit • 
The regression is repeated on each field for many values 
of Tiimit while also varying the velocity range over which 
the data are fit. Not surprisingly the velocity compo¬ 
nents in which the spectra have the largest range in Tb 
yield the most robust fits. 

We find that a and b vary slightly as a function of 
Tiimit and velocity component and that in general the de¬ 
rived slopes and offsets are anti-correlated. With the as¬ 
sumption that errors from removal of stray radiation and 
baseline will result in both over- and under-corrections 

Values for SPIDER but the flanking flelds are similar. 
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Figure 15. Left: Noise map for EBHIS produced from the rms in emission-free channels in each spectrum. Right: Noise map for the 
GHIGLS data adjusted to the same resolution and same gridding. Top to bottom respectively: NEP, NGPLEB, and DRAGO. EBHIS maps 
have been cropped to match the available GHIGLS coverage. Large scale patterns in the noise reflect the different observing strategies. 
Note that the angular scale and the central value and range of the colorbar are different for each map. In NGPLEB emission in the M81 - 
M82 system is sometimes challenging for both GHIGLS and EBHIS, causing some artifacts within about a degree of (/, 6) = (142°, 41°). 


in both datasets, the most robust slopes would be those 
with the lowest absolute offsets. For each velocity com¬ 
ponent we look for the Tnmit that satisfies this criterion 
and also check for the relative constancy of the slope 
over a neighboring range in Tumit- The scale factors 
derived from the Tb and Wm analyses are consistent 
with one another when there is a suitably large range 
in Tb (or VFhi; e.g., Table 3) and likewise the results 
from the three velocity components in each of the three 
fields examined are also consistent. The dispersion of 
the different estimates is 0.030. We adopt a scale factor 
EBHIS/GHIGLS= 1.009 ± 0.011. Before comparing the 
spectral data in the analyses below, we made the minor 


adjustment of all EBHIS survey data to the GHIGLS 
scale using this scale factor. 

The EBHIS and GHIGLS calibrations are consistent to 
within 1%. Erom Boothroyd et al. (2011) the GHIGLS 
calibration has a formal uncertainty of 1 % and there 
we concluded conservatively that “our calibration does 
not have systematic errors that exceed a few percent.” 
We recall that our previous comparison with the LAB 
survey over all GHIGLS fields resulted in a scale factor 
LAB/GHIGLS = 1.0288T0.0012 (Boothroyd et al. 2011). 
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Figure 16. Top panel: Average GHIGLS spectrum for NEP 
(black), average scaled EBHIS spectrum (red), residual scaled 
EBHIS — GHIGLS (blue). Average GHIGLS spectrum peaks at 
5.7 K. Broken lines show the average spectra of the predicted 
stray radiation that has been removed from the data; dashed for 
GHIGLS (black) and dashed and dotted (red) for EBHIS far and 
near side lobes, respectively. Following panels: Average spectrum 
for NGPLEB peaking at 8.1 K and for DRAGO peaking (in the 
IVG) at 3.8 K. 

E.3. Mean and Individual Spectra 

As shown by the average spectra for EBHIS and (con¬ 
volved and regridded) GHIGLS in each field in Figure 16, 
and the difference spectrum, there is remarkable agree¬ 
ment between the two surveys. Also shown are the av¬ 
erage spectra of the predicted stray radiation that have 
been subtracted in producing these spectra. For EBHIS, 
the stray radiation spectrum is calculated in two parts, 
for the near and far side lobes (Kalberla et al. 2010)). 
The GBT has an unblocked aperture, but there is a sig¬ 
nificant spillover sidelobe from the secondary reflector 
(Boothroyd et al. 2011). Obviously, the stray radiation is 
substantial in fields such as surveyed by GHIGLS so that 
without its removal the spectra would be quite different 
and there would be no agreement between the surveys. 
The spectral extent of the stray radiation is larger for the 
far sidelobes where Doppler effects are more significant, 
well into HVG velocities in the case of EBHIS. 

Individual spectra, though noisier, show similar fea¬ 
tures. Figure 17 compares an EBHIS spectrum with the 
corresponding (convolved and regridded) GHIGLS spec¬ 
trum for three distinct lines of sight. The NEP spec¬ 
trum has strong relative emission in the LVG range. 
The NGPLEB spectrum was selected for its relatively 
large GHIGLS predicted stray radiation spectrum. The 
DRAGO spectrum has strong relative emission in the 
IVG range. The difference spectra confirm that the over¬ 
all agreement is very good. 

E.4. Comparisons of Maps of Whi 

In this section we look for evidence in the GHIGLS 
spectra for errors that might be attributable to uncer¬ 
tainties in the baselines that were fit and subtracted or 
to imperfect predictions of the stray radiation spectrum 


Figure 17. Same as Figure 16 but for individual spectra. Top 
panel: NEP spectrum at (99?7, 27?0) peaking at 18.7 K. Following 
panels: NGPLEB spectrum at (137?6,38?1) peaking at 3.2 K and 
DRAGO spectrum at (90?0, 38?7) peaking (in the IVG) at 17.7 K. 

that has been subtracted. The EBHIS survey provides 
an independent basis for this assessment. 

Each of these two contributions is somewhat correlated 
from one channel to the next and thus any errors are 
more easily seen in the line integral VEhl Additionally, 
the velocity dependence of these components suggests 
that any Wm comparison should be made over restricted 
ranges of velocity (e.g., using LVG, IVG, and HVG com¬ 
ponents as defined in Table 2). Gonsequently we produce 
a series of Wm maps corresponding to these velocity com¬ 
ponents. 

In the analysis below of these maps we find some 
evidence for errors relating to the subtraction of base¬ 
lines and stray radiation. These errors are at a low 
level, consistent with the slight increase of cfnui over the 
value expected from noise in the line emission alone (Ap¬ 
pendix D). Overall the GHIGLS and EBHIS maps are in 
good agreement. 

E.4.1. Data 

Figure 18 shows the data available for this investiga¬ 
tion. The field is NEP. The columns (left to right) are for 
integrals over the LVG, IVG, and HVG velocity intervals. 
Unless otherwise indicated, in all of the panels the range 
for the colorbar extends from the minimum to maximum 
of the data values and is usually quite different from one 
panel to another. 

The first row contains the IUhi contributions from the 
baselines that have been fit and subtracted in producing 
the GHIGLS cubes. The second row shows the contri¬ 
butions from the predicted stray radiation for GHIGLS 
that has also been subtracted. In the third row are the 
Hhi maps from GHIGLS already seen in Figure 5. The 
range of the colorbar starts at zero and extends to the 
99 % percentile. For the residual maps in the fourth row 
we have first computed the Hhi maps from the EBHIS 
spectra, adjusted as above to the GHIGLS scale, and 
then subtracted the GHIGLS Hhi maps. Here the color- 
bar is centered on the mean and has a total range that is 
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Figure 18. VUHi-related maps of the NEP field (see text in Section E.4.1 for details). Columns (left to right) are for integrals over the 
LVC, IVC, and HVC velocity intervals. Rows (top to bottom) are from the following data: GHIGLS baselines (subtracted in producing 
GHIGLS Whi), GHIGLS stray radiation (subtracted), GHIGLS Wm (as in Figure 5), VUhi residual from scaled EBHIS - GHIGLS, EBHIS 
far sidelobe (subtracted in producing EBHIS VUhi), and EBHIS near sidelobe stray radiation (subtracted). On the residual map are overlaid 
contours from the GHIGLS stray radiation map (from thinnest to thickest at 25%, 50%, and 75% of the maximum value) for reference 
and orientation. 
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Figure 19. Same as Figure 18 for NCPLEB. In the GHIGLS-related maps, pixels are excluded where it was not possible to fit a simple 
baseline because of emission from other galaxies present in the (frequency-switched) spectrum. Pixels near M81 and M82 where it is not 
possible to fit a baseline are masked and not used in subsequent analysis. A useful indicator of scale is the pattern of the EBHIS multibeam 
system that appears imprinted at one position in the maps of the far sidelobe emission in row five. 
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10 % (20 % for HVC) of that in the Wm map above it. A 
different color table has been used to emphasize the pos¬ 
itive and negative excursions. The fifth and sixth rows 
show the contributions from the predicted stray radiation 
for the far and near sidelobes that have been subtracted 
in producing the EBHIS cubes. The near sidelobe stray 
radiation map is like a blurred version of the original IThi 
map in row three and is therefore roughly correlated with 
the signal. Note that the dynamic range spanned by the 
colorbar is set to be the same as for the far sidelobe map, 
revealing that the corrections for the near sidelobe emis¬ 
sion are generally much smaller. 

Figure 19 contains the corresponding Wni-related 
maps for NCPLEB. 

Even more so than in Figure 15 there are different GPs 
that are obvious in rows one and two and in rows five and 
six. The GPs are related to the scanning strategies of the 
respective surveys and how the observation blocks were 
organized and scheduled. Given their origins the GPs 
are not correlated with the astronomical Wm signal in 
row three, except as already noted for the EBHIS near 
sidelobe stray radiation in row six. 

A strategy for identifying sources of systematic ef¬ 
fects/errors is to search for corresponding telltale pat¬ 
terns like the GPs in the residual maps, i.e., is there 
a “smoking gun” revealing anything amiss with any of 
these corrections? 


E.4.2. Residual Map 

As can be seen from the relative dynamic ranges of 
the maps, the baseline and stray radiation corrections 
that have been subtracted from the GHIGLS and EBHIS 
data are significant. To the extent that they have been 
accurately calculated - and because the data from the 
two telescopes have been set to a common calibration 
scale - the residual maps should appear structureless and 
centred on zero, reflecting solely the combined Hhi noise. 
This is not quite the case. 

The residual map is centered on zero for the LVC com¬ 
ponent in NEP and the I VC and HVC components in 
NCPLEB, and there are broad areas in other maps where 
this is the case too. But when the average residual spec¬ 
trum in Figure 16 is systematically non-zero over a con¬ 
siderable velocity range, for example the negative resid¬ 
ual in the IVC and HVC components in NEP, then this 
will also be the case for the mean of the HAh residual 
map. 

On the residual map we show contours from the 
GHIGLS stray radiation maps to highlight the patchwork 
of largely rectangular patterns relating to the observing 
blocks (the GPs) for reference and orientation. In size 
and orientation the GPs for EBHIS are quite different 
from those for GHIGLS (see also Figure 15). 

The residual map is defined as scaled EBHIS — 
GHIGLS. GHIGLS has been produced from the origi¬ 
nal data minus the corrections in the top two rows and 
so if a correction there were too large it would produce 
an excess in the residual map. Similarly EBHIS has been 
produced from the original data minus the corrections in 
the bottom two rows and so if a correction there were 
too large it would produce a deficit in the residual map. 


E.4.3. Noise in the Line Emission 

The noise in the line emission of the two surveys, of 
order def x [1.0 + Tb(T)/(20 K)] (see Appendix D.3), ac¬ 
cumulates as yNi^Av to contribute fluctuations SW ^ 
0.9 K kms“^ in the residual maps for NEP (0.6 to 
1.3 K kms“^ depending on the field and velocity com¬ 
ponent, the larger values in the more extensive HVC in¬ 
tervals).^^ These fluctuations are small compared to the 
range shown in the residual maps, but do contribute to 
the graininess everywhere. Note that the noise in the line 
emission would not produce any systematic offset in the 
maps. 

E.4.4. Errors Related to Eitted Baselines 

The GHIGLS baselines tend to be fairly stable within 
a given scheduled observing block, and this is reflected 
in the rectangular GPs in the first row of maps. The 
measured standard deviation within a rectangle, typi¬ 
cally 0.5K kms“^ in NEP (0.3 to 1.1 K kms“^ overall), 
contributes to the graininess in the residual maps. The 
measured standard deviation in large smooth areas of 
the residual maps is typically 1.5 K kms“^ in NEP (1.5 
to 1.9K kms“^ overall). Our interpretation is that this 
can be accounted for by the noise in the line emission 
and the baseline fluctuations, including those of EBHIS 
which we have not attempted to quantify. However, we 
note that there are low-amplitude striped patterns in the 
EBHIS and residual maps that can be discerned along the 
distinctively oriented scan lines of EBHIS, with a spacing 
related to the multibeam system, for example in the IVC 
component in NEP and the LVC and IVC components in 
NCPLEB. The measured peak-to-peak amplitude of the 
pattern is small, of order 2K kms“^. 

The GHIGLS baselines do vary from observing block 
to observing block. Just as the GHIGLS baselines are 
fairly stable within the rectangular GPs, any errors in 
the mean offsets within these rectangles are likely to be 
correlated, which could produce a corresponding pattern 
in the residual map. As discussed below, the magni¬ 
tude of such rectangle to rectangle systematic errors is of 
the same order as the above-mentioned fluctuations and 
therefore not normally readily discerned in the residual 
maps given the range of values therein. However, the 
stray radiation correction and the range in the residual 
maps generally decrease steadily from LVC to HVC, mak¬ 
ing HVC IThi maps the most favourable for investigating 
the baseline errors. 

A special case to examine is the large 60 deg^ region in 
NCPLEB ranging in I from 125° and 145° and in b from 
34° to 40° in which there is virtually no HVC emission 
detectable (see the column of HVC maps in Figure 19). 
In the GHIGLS IThi map there are clear GPs relating 
to the observation blocks used, whereas the EBHIS Wm 
map (not shown) is featureless, and so this rectangular 
pattern appears (in reverse) in the residual map. Within 
rectangles in the GHIGLS IThi map the standard devi¬ 
ation is 0.5 K kms“^ whereas the standard deviation for 
the region as a whole is 1.0 K kms“^, from which we de¬ 
duce that the typical dispersion in the mean amplitudes 
of rectangles is 0.9 K km s“^. We also note that even with 

In the GHIGLS Hhi maps themselves the fluctuations would 
be ~ 0.3 K kms“^ for NEP (0.2 to 0.7K kms“^ overall). 
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these clearly visible GPs in the GHIGLS I^hi map, the 
standard deviation is still slightly smaller that it is in the 
smooth EBHIS VEri map (1.1 K hms-^). The GHIGLS 
stray radiation correction map has the same GPs, but 
its overall low level (0.7 ± 0.3K kms“^) and the fact 
that it is uncorrelated with the observed pattern in Wm 
make it an unlikely source. On the other hand there is a 
weak positive correlation of the pattern in IThi with that 
in the baseline correction map. Thus we conclude that 
small observing block to observing block errors in the 
baseline correction with dispersion of about IK kms“^ 
are the likely source of the pattern; this would correspond 
to a rectangle to rectangle dispersion in mean offsets in 
Tb of 8mK over this entire HVG velocity interval. The 
origin of this is unknown but could be related to the gen¬ 
eral problem of interpolating a low-order function over a 
large range. Even when the function is well-constrained 
by the data in many emission-free end channels, the large 
range over which the function is interpolated introduces 
an error. Eor completeness, we note that the average 
VEri in this HVG region is 0.04 K kms-^ for GHIGLS 
and 0.34 K kms“^ for EBHIS. This small difference of 
0.3 K kms“^, the mean of the residual map in this re¬ 
gion, is in accord with the difference in the average spec¬ 
trum (like Eigure 16); this difference could be due to 
small errors in baselines and/or in the stray radiation 
corrections. 

Where there is faint but significant HVG emission de¬ 
tected by GHIGLS in large regions of NEP above 27°, 
and also in DRACO, we find a negative residual, i.e., 
EBHIS underestimates this emission by a few K kms“^. 
We interpret this as a result of the different sensitivi¬ 
ties of the two surveys combined with modelling of the 
baselines. The GHIGLS residual baseline is fit iteratively 
with a third-degree polynomial (Boothroyd et al. 2011), 
whereas the EBHIS baseline is determined using an iter¬ 
ative Gaussian smoothing technique (Winkel et al. 2010) 
applied over their entire bandwidth of 100 MHz. Note 
that the EBHIS data were observed with a frequency¬ 
switching of 3 MHz, but due to the gain varying as a func¬ 
tion of frequency, the data were not reduced using the 
frequency-switching technique (Winkel et al. 2010); thus 
the need for a more sophisticated algorithm to remove 
the baseline. Empirically, the sign of the offset indicates 
that the lower sensitivity (higher noise) of the EBHIS 
spectra may result in a fitted baseline that eliminates 
part of what in GHIGLS is detected as HVG emission 
rather than elevated baseline. Likewise, the GHIGLS 
data baseline fitting could be removing signal that with 
more sensitive observations would lead to the detection 
of resolved structure. 

As mentioned, baseline fitting procedures are not as 
successful at determining and removing a baseline when 
the emission from galaxies dramatically reduces the num¬ 
ber of emission-free channels and/or significantly domi¬ 
nates the shaping of the baseline. 

E.4.5. Stray Radiation 

In Eigure 16 at HVG velocities the average residual 
spectrum contains only low frequency oscillations, con¬ 
sistent with a low-order polynomial used to model base¬ 
lines. This can be contrasted with the sharper oscilla¬ 
tions in the LVG and IVG ranges, which must derive 
from errors other than the baseline. This sharper struc¬ 


ture, also discernible in the residuals of individual spectra 
in Eigure 17, will affect the LVG and IVG Hri resid¬ 
ual maps. The residuals for these components are larger 
than for the HVG. Eurthermore, because of the larger 
dynamic range of the stray radiation corrections in the 
LVG and IVG velocity intervals, these intervals are the 
most favorable for looking for any errors arising because 
of uncertainties in the stray radiation corrections. 

Eor GHIGLS there is no direct spatial correlation be¬ 
tween the stray radiation and Hri because of the un¬ 
blocked geometry of the GBT design. Instead, for in¬ 
termediate latitude fields such as these the amount of 
stray radiation is strongly influenced by the time at which 
any given observation is made because of how the offset 
spillover sidelobe beam pattern (Boothroyd et al. 2011) 
overlaps (or not) with the stronger emission near the 
Galactic plane. This results in a rectangular pattern 
in the stray radiation map relating to the observational 
blocks used. The stray radiation corrections for EBHIS 
are also significant and the maps have quite different 
morphological structure. 

There are some smoking guns. In the NEP residual 
map for IVG in Eigure 18 there are triangular shapes that 
are clearly anti-correlated in amplitude with the EBHIS 
far sidelobe stray radiation correction. To a lesser extent 
there is a positive correlation discernible for LVG. In the 
NCPLEB residual map for LVG in Eigure 19 there are 
rectangular shapes with a clear positive correlation in 
amplitude with the GHIGLS stray radiation correction. 
This persists somewhat in the IVG. 

Attempting some quantification and motivated by this 
suggestion of multiplicative factors, we sought to reduce 
the standard deviation of the residual maps by remov¬ 
ing correlations with various combinations of the stray 
radiation maps. Eor example, in the case of IVG in 
NEP, a factor —0.4 of the EBHIS far sidelobe correc¬ 
tion is indicated (in the case of LVG the factor is -1-0.2). 
However, while this lowers the standard deviation of the 
residual map from 4.6 to 3.7 K kms“^ it is clearly not 
the whole story. There is still a rectangular pattern in 
the revised residual map that reveals the GHIGLS obser¬ 
vational scheduling blocks, but it is not multiplicatively 
related at a significant level (> — O.I). Telltale signs of 
the GPs of both GHIGLS and EBHIS stray radiation cor¬ 
rections in the revised residual maps indicate clearly that 
there are substantial additive and/or subtractive errors 
too, which make any multiplicative factor uncertain. Ac¬ 
counting for these errors is even more important than the 
multiplicative changes to reduce the standard deviation 
in the residual to a level about IK kms“^ that could 
be expected from baseline fitting errors and line noise 
alone. If we suppose that the entire 4.6 K kms“^ dis¬ 
persion is roughly equally apportioned to GHIGLS and 
EBHIS (thus 3.2K kms“^ each), then this amounts to 
about a 10 % fraction of the mean stray radiation correc¬ 
tions (an even small fraction in the case of LVG). How¬ 
ever, we feel that it is probably better to think of the 
errors in absolute terms rather than as a fractional error. 

Eor our NCPLEB LVG example, a factor +0.2 of the 
GHIGLS stray radiation correction is indicated (the stray 
radiation correction has been overestimated). In this 
case this appears to be much larger than the estimated 
VEri stray radiation uncertainty of a factor 0.07 found 
by Boothroyd et al. (2011). But again the standard 
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deviation of the residual map is only lowered from 4.5 
to 3.6 K kms“^ and there are still telltale signs of the 
schedule-related GPs of both GHIGLS and EBHIS stray 
radiation corrections in the original and revised residual 
maps. 

We have carried out a similar comparative study in 
DRACO. In the GHIGLS observations of this field, which 
is a relatively faint extension off the upper right of the 
much larger NEP field, the stray radiation corrections are 
fortuitously at the low end of the values encountered in 
the NEP observations. However, this makes it more diffi¬ 
cult to spot any errors in these corrections. On the other 
hand the EBHIS far sidelobe stray radiation correction is 
quite large and its pattern is discernible particularly on 
the I VC residual map. But because the standard devia¬ 
tion of the residual map is already so low (2.6 K kms“^) 
it is difficult to draw any definitive conclusions, except 
of course that the corrections applied in both surveys 
are very good. This can be appreciated as well in the 
agreement of the mean spectra in Eigure 16. 

F. COMPARISON WITH DATA FROM GASS 

Here we compare our data with that from the southern 
hemisphere-based H I Galactic All-Sky Survey (GASS) 
(McClure-Griffiths et al. 2009; Kalberla et al. 2010; 
Kalberla & Hand 2015) available in the GASS HI data 
release repository. Because of the latitude coverage of 
the GASS data and our choice of primarily northern cir¬ 
cumpolar fields for the surveys with the GBT ACS, we 
have only one field available for this comparison, MC. In 
extracting the GASS cube we adopted an optimal modi¬ 
fied Bessel function mapping of the data onto a 4.'8 grid 
using the gridding parameters recommended by Mangum 
et al. (2007). This is similar but not identical to our 
modified Bessel function mapping of the GHIGLS spec¬ 
tra onto a GLS grid (see Section 3.2). The downloaded 
GASS cube has angular resolution ^ 16' (modified from 
the original average 14.'4 due to our gridding choice) on 
a Cartesian Equatorial grid, 0.82 kms“^ channels, and 
an empirically-determined noise of def ^ 47 mK. 

The GHIGLS data cube was convolved to this lower 
resolution and regridded in all dimensions. This reduces 
the noise from the original 83 mK (Table 1) to an effective 
noise 17 mK. 

Eollowing the same analysis as in Appendix E.2, these 
GASS data appear to be brighter than GHIGLS by a fac¬ 
tor GASS/GHIGLS = 1.062 ±0.019. This is in agreement 
with the conclusion arrived at independently by Kalberla 
& Hand (2015) who compared GASS to EBHIS and LAB 
(LDS). 

On adjusting the GASS data to the GHIGLS scale we 
again find remarkable agreement as demonstrated in the 
spectra in Eigure 20. Also shown there are the stray radi¬ 
ation corrections that have been applied to the GHIGLS 
and the GASS data. 

As in Appendix E.4 we produced a series of Whi maps 
for each velocity component (Table 2) and present them 
in Eigure 21. The residual maps confirm the overall good 
agreement of the two data sets. 

In the HVC part of the spectrum (the largest range, 
adopting —153.4 kms“^ to —32.9 kms“^ from Table 2) 



V (km s 5 



Figure 20. Top panel: Average GHIGLS spectrum for MG 
(black), average scaled GASS spectrum (red), residual GASS — 
GHIGLS (blue). Dashed line spectra are the average stray ra¬ 
diation predictions that have been removed from the data. The 
MG average spectrum peaks at 3.2 K. Bottom panel: Individual 
MG spectrum at (59?4, —81?5) with relatively strong HVG emis¬ 
sion peaking at 2.5 K (LVG peaks at 3.0 K). Note the different 
range as compared to Figures 16 and 17. 


the residual in the average spectrum (Eigure 20) is neg¬ 
ative and therefore so is the residual map for that com¬ 
ponent, which has a mean —4.8K kms“^. The GHIGLS 
stray radiation correction is small by comparison and 
there is no clear imprint of the morphology of the GASS 
stray radiation correction in the Wm residual map. This 
suggests that the origin of the offset might be in the mod¬ 
elling of the baselines. The GASS instrumental baseline 
was removed by an iterative procedure (Kalberla et al. 
2010; Kalberla & Hand 2015) using either a 9th- or 11th- 
order polynomial over a much larger bandwidth than 
used for GHIGLS for which a third-degree polynomial 
could be used (Boothroyd et al. 2011). With the lower 
sensitivity of the GASS measurements the fitting of the 
baseline between the main emission peaks in the profile 
might include what in GHIGLS is detected as emission, 
thus raising the baseline and lowering the apparent sig¬ 
nal, consistent with the sign of the offset. 

Over the IVC range (the smallest, —32.9 kms“^ 
to —8.8 kms“^) the residual in the average spectrum 
changes sign producing a net effect consistent with the 
near-zero mean (±0.7 K kms“^) in the IFhi residual 
map. The dynamic range in both stray radiation correc¬ 
tions is small and in the residual map there is no telltale 
sign of errors in these corrections. However, the location 
and abruptness of the sign change in the residual spec¬ 
tra is suggestive of small errors in the stray radiation 
corrections, rather than in the baselines. 

Over the LVC range (—8.8 kms“^ to ±48.3 kms“^) 
the positive residual in the average spectrum results in a 
slightly positive mean (±2.3 K kms“^) in the Hhi resid¬ 
ual map. There are again no unequivocal signs of errors 
in the stray radiation corrections that have been applied. 

In the GASS data, including the noise map and the 
stray radiation corrections and propagating to the Wm 
residual maps, there is a slight hatching pattern related 
to the orientation of the GASS scans and the multibeam 
system, but this has no bearing on the other considera¬ 
tions here. 
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Figure 21. Same as Figure 18 for MC and comparison with GASS rather than EBHIS. There is only one stray radiation correction cube 
for GASS. Note the areas masked because of the effects of the galaxies MGG-04-02-003, NGG 0045, and NGG 0024 (from left to right). 
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